THE  PREDICTION  OF  PILOT  OPINION 
RATINGS  USING  OPTIMAL  AND 
SUB-OPTIMAL  PILOT  MODELS 

THESIS 

Craig  R.  Edkins,  Captain,  USAF 
AFIT/GAE/ENY/94M-2 

Approved  for  public  release;  distribution  unlimited 


wriv 

lELECTE 
APR  8  81994, 

B 


AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


DTIC  QUALITY  INSPECTED  3 


Wright-Patterson  Air  Force  Base,  Ohio 


AFIT/GAE/ENY/94M-2 


THE  PREDICTION  OF  PILOT  OPINION 
RATINGS  USING  OPTIMAL  AND 
SUB-OPTIMAL  PILOT  MODELS 

THESIS 

Craig  R.  Edkins,  Captain,  USAF 
AFIT/GAE/ENY/94M-2 

Approved  for  public  release;  distribution  unlimited 


94  4  21  061 


94-12280 


March  1994  Master's  Thesis 

THE  PREDICTION  OF  PILOT  OPINION  RATINGS 
USING  OPTIMAL  AND  SUB-OPTIMAL  PILOT  MODELS 


Craig  R.  Edkins,  Captain,  USAF 


Air  Force  Institute  of  Technology 

Wright-Patterson  AFB  OH  45433-6583  AFIT/GAE/ENY/94M-2 


Air  Force  Flight  Dynamics  Directorate 

WL/FIGC  N/A 

Wright-Patterson  AFB  OH  45433-6503 

1 1  SC i"'? t  fvl  i  •  '•  •  <•  \ 


T ) o * *?r ■;  f“;  ' iT;  ■ . '  7y A i  T&.7  :  .  V* V ’ 

APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED 


;3.  A 1 


This  study  details  the  development  of  a  sub-optimal  pilot  model  that  blended  the  classical  and  optimal 
pilot  model  approaches  in  an  attempt  to  achieve  the  advantages  of  each.  This  model  used  a  numerical 
solution  to  the  linear  quadratic  Gaussian  problem  to  find  the  pilot  gain,  lead,  and  lag  values  that 
minimized  a  performance  index  consisting  of  task  error  and  control  usage.  This  development  was 
conducted  in  four  phases.  First,  an  optimal  pilot  model  developed  by  Systems  Technology, 
Incorporated,  was  analyzed  in  detail.  This  analysis  included  a  step-by-step  example  problem  to  clarify 
the  model's  logic  and  an  in-depth  sensitivity  analysis  of  the  model's  parameters.  Second,  a  ground  and 
airborne  evaluation  of  human  pilot  response  was  conducted  using  the  Calspan  variable  stability  Lear  II 
aircraft.  Primary  pilot  response  parameters  were  recorded  and  examined  using  statistical  and  Fourier 
transform  analysis  in  an  attempt  to  provide  insight  into  human  pilot  response.  Third,  a  numerical 
solution  to  the  linear  quadratic  Gaussian  control  problem  that  allows  the  compensator  form  to  be 
predetermined  was  derived.  Finally,  the  sub-optimal  pilot  model  was  developed  and  an  analysis  of  the 
model's  parameters  was  conducted. 
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Preface 


The  purpose  of  this  thesis  was  to  develop  a  pilot  model  that  combined  the 
intuitive  nature  of  classical  pilot  models  with  the  optimal  nature  of  modem  pilot  models. 
This  thesis  is  divided  into  six  sections.  First,  general  background  information  is  given. 
Second,  an  optimal  pilot  model,  developed  by  Systems  Technology,  Incorporated,  is 
analyzed  in  detail.  Third,  a  flight  test  designed  to  provide  insight  into  human  pilot 
behavior  during  ground  and  airborne  tracking  tasks  is  described  in  detail.  Fourth,  a 
numerical  solution  to  the  linear  quadratic  Gaussian  (LQG)  problem  is  derived.  This 
technique  allows  the  compensator  form  to  be  predetermined.  Fifth,  a  sub-optimal  pilot 
model  is  presented.  This  model  uses  the  numerical  LQG  approach  to  restrict  the  optimal 
pilot  model  solution  to  the  classical  pilot  model  form.  Finally,  the  conclusions  and 
recommendations  of  this  thesis  are  summarized. 

This  research  was  accomplished  under  the  joint  Air  Force  Institute  of  Technology 
—  USAF  Test  Pilot  School  program.  I  am  grateful  for  the  unique  opportunity  this 
program  provided.  This  effort  would  not  have  been  possible  without  the  assistance  of 
several  people.  I  would  like  to  thank  my  advisors,  Dr.  Brad  Liebst  and  Lt.  Col.  (Dr.) 
Daniel  Gleason,  for  providing  assistance  and  motivation  during  this  extended  program.  I 
also  wish  to  thank  the  members  of  my  test  management  project  team,  Capt.  Benjamin 
Coffey,  Capt.  Darcy  Granley  (CF),  Capt.  John  Kruzinauskas,  Jr.,  and  Capt.  Mary 
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ABSTRACT 

This  study  details  the  development  of  a  sub-optimal  pilot  model  that  blended  the 
classical  and  optimal  pilot  model  approaches  in  an  attempt  to  achieve  the  advantages  of 
each.  This  model  used  a  numerical  solution  to  the  linear  quadratic  Gaussian  problem  to 
find  the  pilot  gain,  lead,  and  lag  values  that  minimized  a  performance  index  consisting  of 
task  error  and  control  usage.  This  development  was  conducted  in  four  phases. 

First,  an  optimal  pilot  model  developed  by  Systems  Technology,  Incorporated, 
was  analyzed  in  detail.  This  analysis  included  a  step-by-step  example  problem  to  clarify 
the  model's  logic  and  an  in-depth  sensitivity  analysis  of  the  model's  parameters.  This 
study  found  that  when  the  proper  choices  were  made  for  these  parameters,  the  model 
accurately  predicted  Cooper-Harper  ratings.  The  pilot  describing  functions  predicted  by 
this  model,  however,  were  not  consistent  with  the  classical  pilot  model  form. 

Second,  a  ground  and  airborne  evaluation  of  human  pilot  response  was  conducted 
using  the  Calspan  variable  stability  Lear  II  aircraft.  Four  different  pitch  and  four 
different  roll  axis  dynamics  cases  were  evaluated  using  three  different  tracking  tasks. 
Primary  pilot  response  parameters  were  recorded  and  examined  using  statistical  and 
Fourier  transform  analysis  in  an  attempt  to  provide  insight  into  human  pilot  response. 
Except  for  the  presence  of  large  amounts  of  pure  phase  lead  at  high  frequency,  the 
frequency  responses  of  the  pilots  were  consistent  with  the  gain,  lead,  and  lag  form  of  the 
classical  pilot  model.  In  all  cases,  the  pilot  applied  compensation  so  that  the  response  of 
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the  combined  pilot-aircraft  system  resembled  and  integrator  near  the  cross-over 
frequency. 

Third,  a  numerical  solution  to  the  linear  quadratic  Gaussian  control  problem  that 
allowed  the  compensator  form  to  be  predetermined  was  derived.  This  method  used  the 
covariance  matrix  to  compute  the  performance  index,  element  by  element,  and  a 
Nelder-Meade  simplex  search  algorithm  to  find  the  coefficients  of  the  compensator  that 
minimized  the  performance  index.  When  the  desired  compensator  form  was  the  same  as 
the  standard  linear  quadratic  Gaussian  solution,  the  two  methods  produced  identical 
results. 

Finally,  the  sub-optimal  pilot  model  was  developed  and  an  analysis  of  the  model's 
parameters  was  conducted.  This  model  was  restricted  to  single  axis  dynamics  due  to  the 
assumptions  necessary  to  numerically  compute  the  performance  index  value.  By  design, 
the  pilot  describing  functions  predicted  by  the  sub-optimal  pilot  model  were  consistent 
with  the  flight  test  results  and  classical  pilot  modeling  theory.  There  was  an  excellent 
correlation  between  the  performance  indices  and  the  actual  Cooper-Harper  ratings,  but 
the  model  lacked  the  maturity  necessary  for  consistent  Cooper-Harper  ratings  prediction. 
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THE  PREDICTION  OF  PILOT  OPINION  RATINGS 
USING  OPTIMAL  AND  SUB-OPTIMAL  PILOT  MODELS 


1.  INTRODUCTION 


Motivation 

Even  before  Charles  Manley  was  twice  catapulted  into  the  Potomac  River  at  the 
controls  of  the  Langley  Aerodrome,  engineers  were  concerned  about  the  handling 
qualities  of  aircraft.  Early  designs  relied  on  instinct  and  trial  and  error,  often  with 
disastrous  results.  Today,  it  is  critical  that  we  know  as  much  about  an  aircraft's  handling 
qualities  as  early  in  the  design  process  as  possible.  This  analysis  must  include  some  sort 
of  pilot-in-the-loop  evaluation.  The  optimal  and  sub-optimal  pilot  models  are  powerful 
tools  for  performing  this  analysis. 

Pilot  model  development  has  not  kept  pace  with  the  advances  in  flight  control 
design.  Modem  day  performance  and  stealth  requirements  often  dictate  the  geometry  of 
the  aircraft,  and  assume  that  a  digital  flight-control  system  can  be  built  to  make  the 
aircraft  exhibit  good  handling  qualities.  For  this  reason  pilot  model  analysis  is  critical  in 
all  phases  of  an  aircraft's  design.  Unfortunately,  the  handling  qualities  of  these  complex 
flight-control  systems  usually  cannot  be  accurately  predicted  by  currently  accepted  pilot 
models. 
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Historically,  every  advance  in  control  theory  has  led  to  an  attempt  to  understand 
human  pilot  behavior  in  terms  of  the  advance.  The  development  of  classical  control 
theory  after  World  War  II  led  to  a  continuing  effort  to  model  the  human  pilot  with 
transfer  functions.  The  development  of  modem  control  theory  and  linear  quadratic 
Gaussian  techniques  in  the  early  1970s  led  to  the  on-going  attempts  to  model  the  human 
pilot  with  an  optimal  regulator,  filter,  and  estimator.  This  thesis  proposes  a  numerical 
approach  to  pilot  modeling  that  was  made  possible  by  an  exponential  growth  in  computer 
speed  and  availability.  This  numerical  approach  blends  the  classical  and  optimal  pilot 
modeling  theories  in  an  effort  to  achieve  the  advantages  of  each. 


Objectives 

The  overall  objective  of  this  project  was  to  develop  a  pilot  model  that  combines 
the  intuitive  nature  of  classical  pilot  modeling  theory  with  the  optimizing  nature  of 
modem  pilot  modeling  theory.  The  specific  objectives  were: 

1.  Study  existing  classical  pilot  models. 

A.  Implement  the  MIL-STD-1797A  pilot  models  on  MATLAB™1. 

B.  Include  important  aspects  of  these  models  into  the  pilot  model 
developed  for  this  thesis. 

2.  Analyze  the  optimal  pilot  model  developed  by  Systems  Technology, 
Incorporated,  (STI)  for  use  with  Program  CC. 

A.  Supplement  the  guidance  provided  by  STI  for  the  use  of  this  model. 

B.  Conduct  a  sensitivity  analysis  of  the  parameters  used  by  this  model. 

C.  Use  insights  from  this  model  to  develop  the  sub-optimal  pilot  model. 


1  MATLAB  is  a  trademark  of  The  MathWorks,  Inc. 
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3.  Record  and  examine  pilot  response  parameters  during  ground  and  airborne 
tracking  tasks  that  can  be  used  as  a  data  base  for  pilot  model  development 
and  validation. 

A.  Gather  Cooper-Harper  ratings  and  pilot  comments  for  a  range  of 
aircraft  dynamics. 

B.  Using  frequency  response  analysis  techniques,  examine  the 
relationship  between  stick  displacement  and  task  error  for  single  axis 
sum-of-sines  tracking  tasks. 

C.  Use  the  insight  into  human  pilot  response  provided  by  this  experiment 
in  the  development  of  the  sub-optimal  pilot  model. 

4.  Derive  a  numerical  solution  to  the  linear  quadratic  Gaussian  control 
problem  that  allows  the  compensator  form  to  be  predetermined. 

5.  Implement  the  sub-optimal  pilot  model  on  MATLAB™  for  single  axis 
analysis. 

A.  Expand  the  numerical  linear  quadratic  Gaussian  solution  in 
Objective  4  to  include  the  output  disturbance  structure  necessary  for 
the  evaluation  of  single  axis  tracking  tasks. 

B.  Restrict  the  optimal  pilot  model  solution  to  the  classical  pilot  model 
form. 

C.  Compare  the  predictions  of  this  model  with  flight  test  frequency 
response  analysis  data. 

D.  Compare  the  predictions  of  this  model  with  those  of  the  STI  optimal 
pilot  model  and  the  classical  pilot  models  accepted  in 
MIL-STD-1797A. 


Thesis  Overview 

The  following  procedures  were  used  to  accomplish  the  objectives  presented 
above.  Results  are  summarized  at  the  end  of  each  chapter. 

1.  This  thesis  conducted  an  analysis  of  the  crossover  pilot  model  (Reference  17) 
and  the  longitudinal  and  lateral  handling  qualities  models  described  in  MIL-STD-1797A 
(Reference  5).  The  results  of  this  analysis  are  reported  in  Chapter  2  of  this  thesis. 
Additionally,  all  of  the  longitudinal  and  lateral  models  accepted  in  MIL-STD-1797A 
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were  implemented  in  a  Handling  Qualities  Toolbox  for  MATLAB™.  This  toolbox  was 
used  for  much  of  the  classical  pilot  model  analysis  conducted  in  this  thesis.  Copies  of  the 
MATLAB™  programs  (.m  files)  and  user's  guide  are  available  from  the  author. 

2.  The  optimal  pilot  model  developed  by  Systems  Technology,  Incorporated,  was 
analyzed  in  detail.  This  effort  is  detailed  in  Chapter  3.  A  model  overview  is  provided 
and  a  simple  example  is  worked  step-by-step  using  the  program's  logic.  The  sensitivity 
of  the  model  results  to  the  required  parameter  choices  is  also  examined  and 
recommended  values  are  presented. 

3.  A  ground  and  airborne  evaluation  of  human  pilot  response  was  conducted 
using  the  Calspan  variable  stability  Lear  II  aircraft.  The  results  of  this  flight  test  are 
presented  in  Chapter  4.  Four  different  pitch  and  four  different  roll  axis  dynamics  were 
evaluated  using  three  different  types  of  compensatory  tracking  tasks.  Primary  pilot 
response  parameters  were  recorded  and  examined  using  Fourier  transform  analysis  in  an 
attempt  to  provide  insight  into  human  pilot  response.  The  flight  test  data  gathered  during 
this  project  are  maintained  at  the  Flight  Dynamics  Directorate  (WL/FIGC), 
Wright-Patterson  AFB,  Ohio,  and  are  available  for  research  purposes.  A  companion 
report  to  this  thesis,  AFFTC-TLR-93-41  (Reference  7)  serves  as  a  detailed  guide  for  this 
flight  test  data  base  and  provides  an  initial  look  a  the  test  results. 

4.  A  numerical  solution  to  the  linear  quadratic  Gaussian  (LQG)  control  problem 
that  allows  the  compensator  form  to  be  predetermined  is  derived  in  Chapter  5.  This 
method  uses  the  covariance  matrix  to  compute  the  performance  index  value,  element  by 
element.  A  Nelder-Meade  simplex  algorithm  then  finds  the  coefficients  of  compensator 
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that  minimize  the  performance  index.  This  method  is  not  only  useful  in  the  optimal  pilot 
model  problem.  It  may  be  beneficial  for  any  situation  when  reduced  order  compensation 
is  desired. 

5.  The  sub-optimal  pilot  model  is  described  in  Chapter  6.  This  model  uses  the 
numerical  LQG  method  to  restrict  the  optimal  pilot  model  solution  to  the  classical  pilot 
model  form.  The  results  of  this  new  pilot  model  were  compared  with  the  flight  test 
frequency  response  analysis  data  as  well  as  the  predictions  of  the  STI  optimal  pilot  model 
and  the  classical  pilot  models  accepted  in  MIL-STD-1797A. 

6.  Finally,  the  conclusions  and  recommendations  of  thesis  are  summarized  in 
Chapter  7. 

Limitations 

This  thesis  had  several  limitations.  First,  only  single  axis,  compensatory  tracking 
tasks  were  evaluated  or  considered.  Valid  multi-axis  analysis  is  not  possible  without  the 
legitimate  single-axis  framework  this  thesis  was  meant  to  advance.  Second,  the  effects  of 
feel  system  characteristics  were  not  considered.  A  thorough  analysis  of  this  important 
factor  was  beyond  the  scope  of  this  thesis.  Finally,  the  pilots  used  in  the  flight-test 
portion  of  this  study  were  not  entirely  linear,  and  they  did  not  employ  the  ratings  scale 
with  perfect  consistency.  While  every  attempt  was  made  to  reduce  pilot  variations,  these 
restrictions  are  impossible  to  avoid  and  should  be  considered  when  studying  any  handling 
qualities  experiment. 
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2.  Background 


General 

This  chapter  provides  the  background  information  essential  to  this  thesis.  It 
includes  five  sections.  The  first  section,  Handling  Qualities,  discusses  the  Cooper- 
Harper  rating  scale  and  defines  some  important  handling  qualities  terms.  The  next  three 
sections,  The  Pilot  as  a  Linear  Element,  Pilot-in-the-Loop  Analysis,  and  Pilot  Models, 
provide  a  brief  description  of  linear  pilot  modeling  fundamentals.  The  final  section,  Feel 
System  Considerations,  discusses  the  relationship  between  force  and  displacement 
controllers,  and  the  impact  of  each  on  handling  qualities. 

Handling  Qualities 

The  goal  of  this  thesis  is  to  develop  a  pilot  model  that  will,  under  certain 
conditions,  predict  an  aircraft's  handling  qualities.  Cooper  and  Harper  defined  handling 
qualities  as,  " . . .  those  qualities  or  characteristics  of  an  aircraft  that  govern  the  ease  and 
precision  with  which  a  pilot  is  able  to  perform  the  tasks  required  (2:2)."  Central  to  this 
definition  are  the  concepts  of  performance,  workload,  and  task. 

Aircraft  handling  qualities  can  not  be  based  on  performance  alone.  A  pilot  can 
achieve  the  same  performance  for  a  wide  range  of  aircraft  characteristics.  As  the  aircraft 
characteristics  degrade,  however,  this  performance  can  be  attained  only  at  the  expense  of 
a  reduction  in  the  pilot's  capacity  to  perform  other  duties.  Thus,  handling  qualities  also 
depend  on  workload,  defined  to  include  both  mental  and  physical  effort  (2:6). 
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Additionally,  the  same  aircraft  may  exhibit  different  handling  qualities  for 
different  tasks.  For  example,  if  the  task  is  overly  simple  to  perform,  the  poor  handling 
qualities  of  an  aircraft  may  not  be  apparent  to  the  pilot.  Cooper  and  Harper  define  task 
as  the  work  assigned  to  the  pilot.  Task  does  not  include  all  aspects  of  the  mission  or 
intended  use  of  the  aircraft  (2:4). 

The  Cooper-Harper  rating  scale,  shown  in  Figure  2-1  on  the  following  page  and 
described  in  Reference  2,  was  the  basis  for  all  of  the  handling  qualities  ratings  used  in 
this  thesis1.  This  scale  requires  the  pilot  to  answer  a  series  of  questions  concerning 
controllability,  performance,  and  workload  to  arrive  at  a  rating.  Each  rating  is  merely  a 
shorthand  for  a  set  of  handling  qualities  characteristics.  Additionally,  assigning 
Cooper-Harper  ratings  is  a  subjective  process  and  variability  in  the  pilot  ratings  is 
unavoidable.  Because  of  this,  predicting  precise  Cooper-Harper  ratings  with 
mathematical  models  is  problematic  at  best2.  The  pilot  models  examined  in  this  thesis 
were  considered  successful  when  they  consistently  predicted  the  ratings  trends. 

The  Pilot  as  a  Linear  Element 

Pilot  behavior  has  eluded  understanding  since  the  earliest  days  of  aviation.  Not 
only  is  the  human  pilot  non-linear  and  adaptive,  but  his  behavior  varies  over  time.  No 
single  existing  model  can  imitate  the  complexities  of  human  behavior.  Instead, 
simplifying  assumptions  must  be  made  and  the  range  of  the  model's  application 

1  The  term  handling  qualities  rating  and  pilot  opinion  rating  are  interchangeable  in  this  thesis.  Both 
refer  to  the  Cooper-Harper  rating. 

2  Cooper  and  Harper  warn  against  attempting  to  treat  pilot  rating  data  with  mathematical  operations  that 
are  rigorously  applicable  only  to  a  linear,  interval  scale. 
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Figure  2-1.  Cooper-Harper  Rating  Scale 


appropriately  limited.  The  most  significant  assumption  made  throughout  this  thesis  is 
that  the  human  pilot  is  a  linear  element. 

The  response  of  a  non-linear  system  can  be  divided  into  two  parts.  The  first  part 
is  linearly  coherent  with  the  input  and  generally  taking  a  describing  function  form.  The 
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second  part,  the  remnant,  is  not  linearly  coherent  with  the  input.  To  accurately  make  the 
linear  modeling  assumption,  this  remnant  must  be  minimized  (14:4). 

McRuer  found  that  when  modeling  the  human  pilot,  the  remnant  is  primarily  due 
to  four  sources  (14:37). 

1.  Non-Linear  Operations 

2.  Non-Steady  Behavior 

3.  Operator  Response  to  Other  Inputs 

4.  Injection  of  Noise 

The  first  two  sources,  non-linear  operations  and  non-steady  behavior,  result  primarily 
from  the  pilot's  ability  to  learn  and  adapt.  To  minimize  their  effects,  both  the  task  and  the 
pilot  must  be  carefully  chosen. 

The  task  must  be  random  appearing  and  within  the  capabilities  of  the  pilot.  This 
prevents  the  pilot  from  learning  the  task  and  adapting  his  behavior  accordingly.  McRuer 
found  that  the  remnant  also  increased  with  task  bandwidth,  but  the  effect  was  minor  for 
bandwidths  between  1.5  and  4.0  radians  per  second  (15:128).  When  the  bandwidth  of 
the  task  approached  8.0  radians  per  second,  however,  the  pilot  ignored  the  high  frequency 
commands  and  adopted  a  more  open-loop,  non-linear  behavior  (17:246-248). 
Additionally,  the  task  must  have  an  appropriate  length.  It  has  to  be  long  enough  to  allow 
the  pilot  time  to  become  actively  involved,  but  not  so  long  that  he  becomes  fatigued. 
McRuer  found  no  evidence  of  excessive  non-linear  behavior  for  task  lengths  between 
twenty  seconds  and  four  minutes  (17:240). 

Finally,  the  pilot  must  be  carefully  chosen.  McRuer  found  that  using  highly 
trained  and  motivated  subjects  reduced  pilot  variability  and  non-linear  behavior  (15:12). 
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Smith,  in  Reference  24  disagrees.  He  suggests  that  highly  experienced  pilots  use  their 
skill  and  experience  to  avoid  closed  loop  control  when  necessary  (24:20).  While  Smith 
may  be  correct,  it  is  the  goal  of  this  thesis  to  develop  a  pilot  model  that  will  predict  an 
aircraft's  handling  qualities.  The  traditional  evaluator  of  an  aircraft's  handling  qualities  is 
the  experienced  test  pilot. 

The  final  two  remnant  sources,  operator  response  to  other  inputs  and  injection  of 
noise,  are  primarily  affected  by  experimental  procedures.  To  minimize  the  remnant,  the 
task  must  be  clearly  and  carefully  defined.  For  example,  if  a  pure  single  axis  roll 
tracking  task  is  desired,  the  lateral  and  directional  axes  should  be  uncoupled,  or  the  pilot 
will  have  to  divide  his  attention  between  the  two.  These  final  two  remnant  sources  are 
often  mutually  exclusive.  If  a  vision  restriction  device  is  used  so  that  the  only  pilot 
visual  cue  is  task  error,  for  example,  his  response  to  other  inputs  may  be  minimized.  The 
injection  of  noise  may  greatly  increase,  however,  due  to  the  resulting  disorientation. 

The  only  task  examined  in  this  thesis  was  the  compensatory  tracking  task.  For 
this  task,  the  error  between  a  command  and  the  actual  aircraft  condition  was  displayed  to 
the  pilot  and  he  acted  to  minimize  this  error.  This  task  was  considered  representative  of 
several  important  flight  phases,  such  as  air-to-air  and  air-to-ground  gunnery,  aerial 
refueling,  and  approach  and  landing  in  gusty  winds.  The  tasks  used  during  both  the 
analysis  and  the  actual  experiment  in  this  thesis  were  random  appearing  with  a  bandwidth 
of  2.0  radians  per  second.  Highly  trained  pilots  from  the  USAF  Test  Pilot  School  were 
used  as  subjects  in  the  experimental  portion  of  this  thesis.  Additionally,  uncoupled, 
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linear  vehicle  dynamics  were  simulated,  and  the  effects  of  the  dutch  roll  and  phugoid 
modes  of  motion  were  minimized. 


Pilot-in-the-Loop  Analysis 

For  a  large  class  of  tracking  tasks,  the  pilot  acts  to  minimize  the  error  between 
some  aircraft  parameter  and  a  task  command.  Using  pitch  (0)  and  pitch  command  (0C), 
for  example,  this  compensatory  tracking  task  is  illustrated  in  Figure  2-2. 


0c 


*  Y_ 

8s  or  Fs 

v 

*  YP 

0 


Figure  2-2.  Compensatory  Tracking  Task  Block  Diagram 


where 

0  =  Pitch  Angle  Yp  =  Pilot 

0C  =  Pitch  Command  Yc  =  Aircraft 

e  =  Task  Error  Fs  =  Stick  Force 

8,  =  Stick  Deflection 

Through  block  diagram  manipulation  this  can  be  re-drawn  such  that  the  pilot  senses  the 
task  error  and  applies  compensation  to  minimize  it  as  shown  in  Figure  2-3.  This  form 
will  be  used  in  the  remainder  of  this  thesis  for  all  model  development.  The  salient 
question  is,  "  What  is  Yp  ?” 
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Figure  2-3.  Compensatory  Tracking  Task  Block  Diagram  -  Feedback  Compensation 


In  1965  McRuer  and  Krendel  conducted  a  controlled  measurement  of  human  pilot 
behavior  using  frequency  response  methods.  Using  a  pitch  axis  tracking  task  generated 
on  a  laboratory  oscilloscope,  they  found  that  for  very  simple  dynamics,  the  quasi-linear1 
pilot  can  be  represented  by  the  following  describing  function  (15: 14-17). 


Yp  =  Kp  •  e-»x 


Gain  Pure 
Delay 


(Tl-Jo  +  I)  [TKj®  +  l] _ 1 _ 

(Trja+l)  [fK  •>  + 1]  (Tm  ,j(a  + 1} .  f  (g)2  +  feo  +  1 

v  ^  v  L 


Series  Very  Low 

Equalization  Frequency 

Lag-Lead 


Neuromuscular 

System 


where 


j©  =  Laplace  Variable  for  Random  Input2 

Yp  =  Pilot  Describing  Function 

Kp  =  Pilot  Gain 

x  =  Pilot  Delay 

Tl  =  Pilot  Lead  Time  Constant 

T,  =  Pilot  Lag  Time  Constant 

Tk  =  Very  Low  Frequency  Pilot  Lead  Time  Constant 

T'k  =  Very  Low  Frequency  Pilot  Lag  Time  Constant 

TN1  =  Neuro-Muscular  Time  Constant 

©„  =  High  Frequency  Neuro-Muscular  Natural  Frequency 

Q,  =  High  Frequency  Neuro-Muscular  Damping  Ratio 


1  McRuer  used  the  term,  quasi-linear,  because  the  pilot  can  be  modeled  as  a  linear  element  only  if 
specific  conditions  are  satisfied. 

2  jw  is  used  instead  of  the  general  Laplace  variable,  s  =  a+jco,  to  emphasize  that  this  equation  is  strictly 
valid  only  in  the  frequency  domain  with  continuous,  random-like  inputs.  It  should  not  be  used  for  system 
responses  to  deterministic  inputs  such  as  step  commands. 
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Limiting  the  frequency  range  of  interest  to  0.1  to  10  radians  per  second,  the  influence  of 
the  very  low  frequency  lag-lead  is  minimal.  Additionally,  the  neuro-muscular  system  can 
be  modeled  as  a  first  order  lag  as  shown  in  Equation  2-2  (16:29). 


_ i _  M  i 

-yeo  + 1)  •  f (&1  +  +  ll  ~  7W“  +  1 


if  7V*7V,+f^  (2-2) 


where 

Tn  =  First  Order  Approximation  of  Neuro-Muscular  Time  Constant 

and  typical  values  are  (15:171): 

TN1  =  0.1  C  =  0.12 

©„  =  16.5  Tn  *  0.115 

Making  these  assumptions  leaves  the  classical  pilot  model  shown  in  Equation  2-3. 


Yp  =  Kp  ■  e-J"' 


_ (77,  j®  +  1) _ 

(77  j®  + 1)  •  (TN  jo  + 1) 


(2-3) 


As  shown  in  this  equation,  the  pilot  can  be  modeled  by  a  gain,  a  lead,  a  lag,  a  delay,  and  a 
first  order  muscular  lag. 

Using  the  classical  pilot  model  and  including  measurement  and  muscular  noise 
the  block  diagram  of  Figure  2-3  can  be  re-drawn  as  in  Figure  2-4  on  the  following  page. 
As  shown  in  this  figure,  the  pilot  makes  an  imperfect,  or  noisy  (£c),  measurement  of  the 
task  error  and  applies  a  delayed  compensation  to  minimize  this  error.  The  desired  stick 
deflection  (8C)  is  corrupted  by  both  the  muscular  lag  and  a  random  muscular  noise  (£,m ). 
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Figure  2-4.  Classical  Pilot  Model  Block  Diagram 


where 

0  =  Pitch  Angle  4<>  =  Observation  Noise 

0C  =  Pitch  Command  4m  ~  Muscular  Noise 

e  =  Task  Error  5C  =  Commanded  Stick  Deflection 

Kp  =  Pilot  Gain  5S  =  Stick  Deflection 

Tl  =  Pilot  Lead  Time  Constant  F,  =  Stick  Force 

Tj  =  Pilot  Lag  Time  Constant  Yc  =  Aircraft  Transfer  Function 

x  =  Pilot  Delay  jco  =  Laplace  Variable 

Tn  =  Muscular  Time  Constant  c/r  =  Output/Input 

While  pilot  lead  makes  practical  sense,  pilot  lag  is  not  intuitively  obvious.  Smith 
proposed  a  more  intuitive  approach.  In  Reference  24  Smith  theorized  that  the  pilot 
measures  both  error  and  error  rate  and  applies  a  gain  to  each,  resulting  in  lead 
compensation.  Thus,  the  larger  the  error  or  error  rate,  the  larger  the  compensating  stick 
movement,  but  only  to  a  point.  As  the  error  or  error  rate  grow  large,  the  pilot  will  no 
longer  proportionally  increase  the  stick  deflection  by  the  same  amount.  In  other  words 
the  pilot  will  apply  washout  to  the  error  and  error  rate  signals.  It  is  this  washout  of  the 
error  and  error  rate  signals  that  produces  the  lag  in  Equation  2-3,  not  some  mentally 
computed  lag  compensation  (24:37-40).  This  can  be  drawn  as  shown  in  Figure  2-5. 
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Figure  2-5.  Multi-Loop  Pilot  Model  Block  Diagram 

where 

T*,,*  Washout  Time  Constant  in  Error  Feedback  Loop 
Tw  -  Washout  Time  Constant  in  Error  Rate  Feedback  Loop 
Kp,  -  Error  Feedback  Gain 
=  Error  Rate  Feedback  Gain 

It  is  interesting  to  note  that  the  multi-loop  model  in  Figure  2-5  and  the  classical 
pilot  model  in  Figure  2-4  are  equivalent  if  the  washout  time  constants  in  the  error  and 
error  rate  feedback  loops  are  equal  (Twp=  Tw).  To  make  this  illustration  simpler,  the 
identical  forward  paths  in  each  figure  were  replaced  with  the  variable,  G.  Using  block 
algebra,  the  closed-loop  transfer  function  for  Figure  2-4  is: 


c  _ _ G  •  (7/  -jd)  +  1) _ 

r  ~  Ti  j(a  + 1  -  G  Kp  •  (Tl  j(o  + 1)  •  T 


(2-4) 


while  the  closed-loop  transfer  function  of  Figure  2-5  is: 

c  _ _ G  •  (T^p  j(o  + 1)  •  (Tht  j<n  + 1) _ 

r  (Tyvp  jm  +  l-G  ■  Kpi  ■  e"^'T)  •  ( I ^  j(o  + 1)  -  G  •  (7^  j(o  +  1)  j(o  ■  Kpl  •  T 
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If  the  washout  time  constants  in  both  the  error  and  error  rate  feedback  loops  are  equal 
(Tw,,=  TWI=  Tw),  then  Equation  2-5  reduces  to: 


c _ G  ■  (Tw  -j( o  +  1) _ 

r~  Twj<a  +  \-G- Kpl  ■  e-*°"  - Gja> ■  Kp2  ■ 


(2-6) 


which  can  be  rewritten  as: 

c _ G  *  (Tw  'j®  + 1) _ 

Twja>  +  l-GKpl(^ja>  +  l) -er*“ 


(2-7) 


This  transfer  function  is  equivalent  to  that  in  Equation  2-4  found  using  the  classical  pilot 
model.  Several  important  insights  can  be  gained  from  comparing  Equations  2-4  and  2-7. 
The  classical  pilot  modeling  approach  and  the  multi-loop  approach  are  identical  if: 

1 .  The  washout  time  constant  in  both  the  error  and  error  rate  feedback 
loops  of  the  multi-loop  model  are  equal  to  the  lag  time  constant  of  the 
classical  pilot  model  (T^  =  Tw=  Tw=  Tj). 

2.  The  gain  on  error  feedback  in  both  models  are  equal  (K^,  =  K^). 

3.  The  ratio  of  error  rate  to  error  feedback  gain  in  the  multi-loop  model 
equals  the  lead  time  constant  in  the  classical  pilot  model 

(K*/  K,.-T,.). 

Pilot  Models 

There  are  currently  three  broad  categories  of  pilot  models.  The  first  category, 
open  loop  pilot  models,  uses  the  aircraft  response  to  open  loop  commands  to  predict 
handling  qualities.  These  models  are  based  on  statistical  fits  to  past  data  and  make  no 
attempt  to  model  human  pilot  behavior.  The  second  category,  classical  pilot  models, 
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model  the  pilot  as  a  gain,  lead,  lag,  and  delay  as  shown  in  Equation  2-3.  The  final 
category,  optimal  pilot  models,  model  the  pilot  so  that  some  performance  index  is 
minimized.  These  models  generally  take  the  form  of  a  linear  quadratic  regulator  and 
Kalman  filter. 

The  majority  of  pilot  models  in  use  are  open  loop  models.  Because  these  models 
do  not  directly  consider  human  pilot  behavior,  they  are  comparatively  simple  to  develop. 
Five  of  the  six  pitch  axis  and  all  of  the  roll  axis  models  in  MIL-STD-1797A  are 
open-loop  models  (5:171-256,  377-427). 

Classical  and  optimal  pilot  models  face  three  obstacles.  First,  a  form  for  the  pilot 
must  be  chosen.  For  example,  the  classical  pilot  models  use  gain,  lead,  lag,  and  delay. 
Second,  the  constants  in  the  pilot  model  forms  must  be  found.  The  optimal  pilot  model 
computes  these  values  so  that  a  performance  index  is  minimized.  Finally,  the  resulting 
pilot  describing  function  must  be  related  to  the  aircraft's  handling  qualities  through  some 
metric. 

The  primary  advantage  of  classical  pilot  models  is  that  the  form  used  by  these 
models  is  based  on  experimental  results.  The  gain,  lead,  lag,  and  delay  behavior  of  the 
human  pilot  during  compensatory  tracking  tasks  is  well  documented  (References  14,  15, 
16, 17,  and  21).  Choosing  values  for  these  variables  and  then  using  them  to  predict  a 
handling  qualities  rating,  however,  can  be  quite  difficult. 

Classical  pilot  model  parameter  selection  does  have  some  experimental  basis. 
McRuer  found  that  the  pilot  uses  gain,  lead,  and  lag  so  that  just  within  and  beyond  the 
task  bandwidth  the  combined  aircraft-pilot  system  has  a  slope  of  -20  dB  per  decade  and 
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adequate  stability  margins.  Changes  in  aircraft  gain  are  offset  by  changes  in  pilot  gain  so 
that  the  crossover  frequency  is  largely  invariant  (17:234-238).  McRuer  also  found  that 
the  generation  of  second  order  lead  was  difficult  or  impossible  with  only  visual  random 
appearing  inputs.  While  several  of  the  systems  he  tested  could  have  used  second  order 
compensation  to  advantage,  this  behavior  was  not  noted  (14:218).  Finally,  McRuer 
measured  pilot  delay  between  0.2  and  0.3  seconds  (14:217). 

The  crossover  pilot  model  used  this  behavior-based  approach  to  form  a  series  of 
rules  for  selecting  pilot  gain,  lead,  and  lag  (Reference  17).  The  problem  of  accurately 
and  consistently  relating  these  variables  to  a  handling  qualities  rating  still  remains.  As  a 
result,  classical  models  are  not  in  wide  use.  The  Neal-Smith  model  is  currently  the  only 
classical  pilot  model  accepted  in  MIL-STD-1797A. 

Optimal  pilot  models  are  based  on  the  assumption  that  the  well-trained, 
well-motivated  human  controller  will  perform  in  a  near  optimal  manner  subject  to  certain 
internal  constraints.  Optimal  pilot  models  have  two  advantages  over  classical  pilot 
models.  First,  selection  of  the  describing  function  constants  is  not  difficult.  The  values 
are  chosen  so  that  the  linear  quadratic  Gaussian  performance  index  is  minimized. 

Second,  relating  the  pilot  describing  function  to  a  handling  qualities  rating  is  more 
straight  forward  than  with  classical  pilot  models.  Usually,  the  performance  index  is 
directly  related  to  a  handling  qualities  rating. 

The  optimal  pilot  model  suffers  from  two  deficiencies.  First,  optimal  pilot 
models  are  complicated,  both  in  the  order  of  the  predicted  pilot  describing  function  and 
in  their  implementation.  Due  to  the  nature  of  linear  quadratic  Gaussian  theory,  the 
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predicted  pilot  describing  function  has  an  order  of  at  least  twice  that  of  the  aircraft 
dynamics.  This  is  not  consistent  with  the  gain,  lead,  lag,  and  delay  behavior  observed  by 
McRuer.  Additionally,  implementing  the  internal  constraints  of  the  human  pilot,  delay 
and  muscular  lag  for  example,  within  a  standard  linear  quadratic  Gaussian  structure  is 
complicated.  The  second  deficiency  is  over  parameterization.  When  implementing  an 
optimal  pilot  model,  the  engineer  has  to  make  numerous  parameter  choices,  such  as  noise 
intensities,  performance  index  weightings,  and  forcing  function  form  and  intensity.  The 
predicted  handling  qualities  rating  and  pilot  describing  function  is  normally  sensitive  to 
these  choices. 

Feel  System  Considerations 

It  is  difficult  to  imagine  any  pilot-in-the-loop  analysis  without  considering  control 
stick  characteristics.  Controller  sensitivity,  or  stick  force  gradient,  can  have  a  huge  effect 
on  an  aircraft's  handling  qualities.  In  one  experiment,  a  Cooper-Harper  rating  was 
changed  from  a  seven  to  a  two  only  by  changing  the  roll -rate  sensitivity  (18:19).  Control 
stick  dynamics  and  control  stick  type  —  force  versus  position  command  --  are  also 
important  factors  in  an  aircraft’s  handling  qualities. 

The  position  command  feel  system  dynamics  include  the  spring,  mass,  and 
damper  characteristics  that  translate  the  pilot's  stick  force  input  into  stick  position.  There 
are  two  different  approaches  to  analyzing  this  type  of  feel  system.  It  can  be  considered 
as  another  flight  control  system  filter,  or  it  can  be  treated  as  a  unique  dynamic  element. 

If  the  feel  system  is  treated  as  another  flight  control  system  filter,  the  stick  dynamics  are 
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included  in  the  analytical  model  as  a  linear  element  between  stick  force  and  aircraft 
response  as  shown  in  Figure  2-6  (18:36). 


Figure  2-6.  Position  Command  Feel  System 

If  the  feel  system  is  treated  as  a  unique  element,  it  is  assumed  that  the  pilot  uses  his  direct 
access  to  the  input  force  and  the  output  stick  displacement  to  form  an  inner  loop  around 
the  feel  system.  Thus,  the  stick  dynamics  become  internal  to  the  pilot  and  do  not  need  to 
be  included  in  the  analytical  model  (23:582-583).  MIL-STD-1797A  advocates  the  latter 
approach  (5:172). 

The  force  command  feel  system  acts  in  parallel  as  shown  in  Figure  2-7  (18:36). 


To  Flight 
Control  System 

Figure  2-7.  Force  Command  Feel  System 

The  traditional  approach  is  to  exclude  the  feel  system  dynamics  for  force  controllers 
since  they  act  in  parallel.  Without  the  feel  system  dynamics,  the  system  will  exhibit  less 
phase  lag  and  therefore  better  closed-loop  performance  (18: 14). 

Mitchell  found  that  including  the  feel  system  as  an  equivalent  delay,  regardless  of 
the  feel  system  type,  improved  the  correlation  between  handling  qualities  ratings  and  feel 
system  dynamics.  In  other  words,  an  aircraft  with  a  "bad"  feel  system  should  receive 
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poor  ratings  whether  position  or  force  command  is  used.  The  systems  could  be  modeled 
as  an  equivalent  delay  due  to  the  neuromuscular  reaction  of  the  pilot  to  both  types  of  feel 
systems.  Mitchell  emphasized  that  the  effect  of  the  feel  system  was  not  identical  to  a 
pure  time  delay.  It  was  possible,  however,  to  estimate  the  effects  of  the  feel  system  on 
pilot  ratings  by  considering  it  as  such  (18:24). 

This  thesis  only  studied  position  command  systems,  and  the  feel  system  dynamics 
were  generally  analyzed  as  a  linear  element  in  the  command  path.  An  in-flight  analysis 
of  both  position  and  force  command  systems  was  not  feasible  because  of  the  numerous 
additional  flight  test  combinations  it  would  have  required. 
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3.  The  Optimal  Pilot  Model 


General 

While  the  gain,  lead,  lag,  and  delay  form  of  the  classical  pilot  model  is  well 
documented,  choosing  values  for  these  variables  and  then  using  them  to  predict  a 
handling  qualities  rating  remains  a  daunting  task.  The  optimal  pilot  model  was 
developed  in  an  attempt  to  overcome  these  deficiencies.  The  optimal  pilot  model  is  based 
on  linear  quadratic  Gaussian  (LQG)  stochastic  control  theory  which  appears  intuitively 
suited  to  the  solution  of  the  human  pilot  modeling  problem.  The  human  pilot  adapts  his 
control  strategy  so  that  the  trade-off  between  performance  and  workload  results  in  the 
optimal  handling  qualities  rating  (6:7).  Likewise,  the  optimal  pilot  model  finds  the 
control  strategy  that  minimize?  a  performance  index  consisting  of  average  tracking  error 
( performance )  and  control  usage  (workload).  This  performance  index  can  be  used  to 
estimate  a  handling  qualities  rating  based  on  statistical  fits  to  past  data. 

This  chapter  focuses  on  an  optimal  pilot  model  developed  by  Systems 
Technology,  Incorporated  (STI),  for  use  with  Program  CC  (Reference  25).  The  STI 
model  was  selected  for  several  reasons.  First,  it  incorporates  nearly  every  important 
aspect  of  other  optimal  pilot  models,  making  it  a  good  candidate  for  study.  Second,  it 
can  be  implemented  on  the  personal  computer  and  is  therefore  widely  available.  Finally, 
this  model  has  had  some  success  in  predicting  Cooper-Harper  ratings,  but  lacks 
parameter  selection  guidance. 
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This  chapter  is  divided  into  three  sections.  The  first  section  gives  a  brief 
overview  of  optimal  pilot  modeling  theory.  The  second  section  analyzes  the  STI  optimal 
pilot  model  in  detail.  This  analysis  includes  a  description  of  the  model,  a  step-by-step 
example  problem,  and  a  detailed  sensitivity  analysis.  Finally,  the  conclusions  and 
recommendations  of  this  chapter  are  summarized. 

Background 

Kleinman,  Baron  and  Levison  made  the  first  integrated  attempt  to  describe  the 
behavior  of  the  human  pilot  in  the  framework  of  optimal  control  theory  in  1970 
(Reference  10).  The  theory  involved  has  not  changed  significantly  since  the  Kleinman 
study  was  completed.  The  basic  assumption  implicit  in  the  optimal  pilot  model  is  that 
the  well-trained  and  motivated  pilot  behaves  in  an  optimal  manner  subject  to  his  inherent 
limitations  (8:464).  These  limitations  include  neuromuscular  lag,  internal  motor  noise, 
perceptual  threshold,  and  processing  delay. 

A  block  diagram  of  the  standard  optimal  pilot  model  is  shown  in  Figure  3-1 
(6:15).  The  aircraft  dynamics  are  represented  by  a  linear  state  vector  and  vector-matrix 
equations.  The  task  forcing  function  is  input  as  a  disturbance  and  modeled  as  white 
Gaussian  noise  shaped  by  coloring  filters.  Because  the  human  pilot  is  able  to  obtain  both 
displacement  and  rate  information  from  a  single  display,  the  feedback  vector  normally 
consists  of  the  error  and  error  rates  of  the  displayed  variables.  If  the  error  or  error  rates 
are  smaller  than  a  threshold  value,  the  pilot  will  neither  observe  them  nor  take  any 
corrective  action. 


3-2 


Pilot's  Control 
Output 


c 

•G  © 

iOuJ1 

is 

Jo  o 
OZ 


-  V 

2.2 

IS 


3-3 


Figure  3-1.  Block  Diagram  of  Optimal  Pilot  Model 


Pilot  induced,  uncorrelated  white  noise  is  added  to  each  observation.  This  noise 
represents  the  nonlinear  remnant  of  the  pilot  as  well  as  his  random  error  in  observing  the 
displayed  variables.  If  more  than  one  variable  and  its  derivative  are  observed,  sampling 
must  also  be  considered  (6: 16). 

The  many  internal  time  delays  inherent  to  the  human  pilot  are  modeled  by  a 
single  perceptual  delay.  A  Kalman  filter  and  least-squares  predictor  are  used  to  produce  a 
non-delayed,  optimal  estimate  of  the  observation  vector.  The  linear  quadratic  regulator 
(LQR)  gains  are  then  determined  to  minimize  an  appropriate  performance  index1.  For 
the  optimal  pilot  model,  this  performance  index  generally  takes  the  form  of  Equation  3-1 . 


J=  J  [xT(t)Qx(t)  +  uT(t)Ru(t)]dt  (3-1) 

o 


where 

x  (t)  =  State  Vector 
J  =  Performance  Index 
Q  =  State  Deviation  Weighting 
Matrix 


u  (t)  =  Control  Rate  Vector 
t  =  Time 

R  =  Control  Rate  Weighting  Matrix 


The  first  part  of  this  integral  (x?Q  x)  is  directly  related  to  performance,  while  the  second 
part  ( iiTRu )  represents  physical  and  mental  workload  (8:475).  The  state  weightings  can 
be  difficult  to  pick  a  priori.  Hess  found  success  by  choosing  the  weighting  coefficients 

as  the  squares  of  the  reciprocals  of  the  maximum  allowable  deviation  of  the  state 

variables  (8:470).  These  weightings  can  also  be  used  to  make  the  model  task  dependent 

(6:11). 


1  Background  information  on  control  rate  weighting  is  provided  in  Appendix  A.  For  more  detailed 
information  on  LQG  theory  consult  References  13  and  20. 
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The  neuromuscular  lag  is  not  normally  explicitly  modeled.  Kleinman,  Baron,  and 
Levison  showed  that  the  neuromuscular  dynamics  did  not  need  to  be  introduced  if  the 
performance  index.  Equation  3-1,  includes  a  penalty  on  control  rate  instead  of  control 
(8:465).  The  neuromuscular  noise,  or  motor  noise,  and  the  observation  noise  have 
spectral  densities  proportional  to  the  mean  squared  operator  output.  Thus,  they  must  be 
found  in  an  iterative  fashion. 

Finally,  the  performance  index  is  used  to  estimate  a  handling  qualities  rating 
based  on  statistical  fits  to  experimental  data  (6:24).  This  relationship  normally  accounts 
for  the  strength  of  the  forcing  function  as  well  as  the  forcing  function  bandwidth. 

STI  Optimal  Pilot  Model 

This  section  presents  a  detailed  analysis  of  the  optimal  pilot  model  developed  by 
Systems  Technology,  Incorporated,  for  use  with  Program  CC.  Reference  25  explains  the 
theory  behind  this  model.  This  section  is  intended  to  supplement  Reference  25  by 
providing  an  extended  discussion  of  the  model  as  well  as  the  user-selected  model 
parameters.  The  discussion  that  follows  is  divided  into  three  parts.  The  first  section 
provides  an  overview  of  the  model.  The  second  section  follows  the  model's 
computational  flow  using  an  integrator  example.  The  final  section  contains  a  detailed 
sensitivity  analysis  of  the  modeling  parameters  and  provides  guidance  for  their  selection. 

Model  Overview.  The  flow  diagram  for  the  STI  optimal  pilot  model  is  shown  in 
Figure  3-2  on  the  following  page  (25:7).  This  model  was  only  designed  to  handle  one 
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Where 

Tn  -  Muscular  Time 
Constant 

g  ->  Control  Rate 
Weighting 

q  -  State  Deviation 
Weighting 

x  -  Pilot  Delay 

Vw  -  Intensity  of 
Driving  Noise 

p*  -  Observation 
Noise  Ratio 

p,  =  Motor  Noise 
Ratio 

T,  =  Visual  Indiffer¬ 
ence  Threshold 

f  =  Fractional  Atten¬ 
tion  Parameter 

Vj,  =  Observation  Noise 
Intensity 

V,  =  Motor  Noise 
Intensity 

=  Standard  Devia¬ 
tion  of  Error  or 
Error  Rate 

=  Standard  Devia¬ 
tion  of  Control 
Rate 

E  -  Error  Function 

J  “  Performance 
Index 

Yr  =  Predicted  Pilot 
Transfer  Function 


Figure  3-2.  Computational  Flow  for  the  STI  Optimal  Pilot  Model 
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feedback  channel  and  its  derivative.  If  there  is  more  than  one  channel,  each  must  be 
analyzed  separately.  This  model  can  be  divided  into  three  parts. 

First  the  aircraft  dynamics  are  input  in  either  transfer-function  or  state-space  form 
and  appended  with  the  necessary  coloring  filters1.  Because  this  model  allows  the 
feedback  of  only  one  variable  and  its  derivative,  single-input-single-output  (SISO) 
aircraft  and  shaping  filter  dynamics  must  be  used.  Once  these  dynamics  are  entered,  the 
aircraft  states  are  augmented  with  the  driving  noise  filter  states  in  one  of  three  different 
manners.  The  noise  can  be  injected  at  the  aircraft  output,  at  the  aircraft  input,  or  in  the 
middle  of  the  aircraft.  Finally,  the  state-space  representation  of  the  augmented  aircraft  is 
configured  so  that  error  and  error  rate  are  used  for  feedback. 

Second,  the  linear  quadratic  regulator  problem  is  solved  so  that  the  neuromuscular 
dynamics  are  modeled.  The  user  enters  the  state  deviation  penalty,  q.  This  weighting 
applies  only  to  the  error  channel  since  the  penalty  on  error  rate  is  always  zero  in  this 
model.  Integrators  are  then  appended  to  the  aircraft  dynamics.  Finally,  the  weighting  on 
control  rate,  g,  is  iterated  until  the  neuromuscular  dynamics  are  modeled.  Because  the 
STI  optimal  pilot  model  is  restricted  to  single-input-single-output  dynamics,  q  and  g  are 
scalar  weightings  and  are  used  in  the  STI  documentation  instead  of  the  standard 
weighting  matrices,  Q  and  R,  in  Equation  3-1. 


1  The  tracking  task  forcing  function  must  be  modeled  as  a  white  Gaussian  noise  shaped  by  coloring 
filters. 
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Finally,  the  Kalman  filter  problem  is  solved.  The  user  must  enter  several 
parameters: 

t  Pilot  Delay 

Vw  Intensity  of  the  Driving  Noise 

Pyi »  Py2  Observation  Noise  Ratios  for  Error  and  Error  Rate 

pu  Motor  Noise  Ratio 

Vy,j ,  Vy2i  Intensity  of  the  Observation  Noise  —  Starting  Values 
Vui  Intensity  of  the  Motor  Noise  --  Starting  Value 

TyI ,  Ty2  Visual  Indifference  Thresholds 

/  Fractional  Attention  Parameter 

The  starting  values  for  the  noise  intensities  are  used  to  find  the  Kalman  filter  solution. 
From  this  solution  ayi  and  ou ,  the  standard  deviations  of  y;  and  u,  are  determined.  The 
no  >c  intensities  are  then  iterated  until  the  equations  for  Vyi  and  Vu  shown  in  Figure  3-2 
are  satisfied.  For  multi-axis  tasks  this  entire  model  must  be  repeated  and  the  fractional 
attention  parameter  manually  iterated  until  the  combination  that  produces  the  minimum 
performance  index  is  found.  Finally,  the  model  computes  the  performance  index  value 
and  a  pilot  describing  function.  The  performance  index  can  be  used  in  Equation  3-2  to 
estimate  a  Cooper-Harper  rating  (19:40). 

Rating  =  5.5 +  3.7  •  log10  —J  (3-2) 

\ai  -GiiJ 

where 

oc  =  Forcing  Function  Root  Mean  Square  Error  Amplitude 
(Dw  =  Forcing  Function  Bandwidth 

Integrator  Example.  The  easiest  way  to  understand  how  this  model  works  is  to 
follow  the  computational  flow  using  a  simple  example.  This  example  will  also  give 
insight  into  the  parameter  selection  process  that  will  be  discussed  in  the  next  section. 
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The  simplest  choice  is  to  assume  the  aircraft  dynamics  can  be  modeled  by  an  integrator 
and  the  forcing-function  can  be  modeled  by  a  first  order  lag.  The  following  values 
parallel  a  sample  problem  presented  in  Reference  25. 


State  Space  Realization 
Yc  =  ^  =  1  =>  [*.]  =  [0]xi+[l]8e 


y  =  6  =  [ljxi  (3-3) 

Yw  =  it =  7TT  =*  l*2]  =  [~2]X2  +  [11  W” 

y=Wc  =  [j 8T]x2  (3-4) 

where 

Yc  =  Aircraft  Dynamics 

Yw  =  Tracking  Task  Forcing  Function 

x,  =  Pitch  State  (0) 

x2  =  Coloring  Filter  State 

Ww  =  Zero  Mean  White  Gaussian  Noise 

Wc  =  Colored  Noise 


The  filter  states  are  now  augmented  to  the  aircraft  states  in  one  of  three  ways. 

The  advantages  and  disadvantages  of  each  method  will  be  discussed  in  the  sensitivity 
analysis  section  of  this  chapter.  The  first  option  is  to  inject  the  noise  at  the  aircraft  output 
according  to  Figure  3-3. 


Figure  3-3.  Noise  Injected  at  Aircraft  Output 
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Because  the  coloring  filter  and  aircraft  dynamics  operate  in  parallel,  their  states  can 
simply  be  appended  together  as  shown  in  Equation  3-5. 


*1 

’  0 
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Xl 

4. 

’  1 

0  " 

_  *2  . 

0 

-2 

.  *2  . 

T 

0 

1 

6« 

Ww 


(3-5) 


The  output  vector,  y,  contains  the  error  and  error  rate  where  this  error  can  be  written  as: 


e=  Wc  +0 


e=Wc  +  6 


but 


0  =  X\ 

6  =ii  =  8, 


and 


wc  -  JfTsx2 


Wc  =  J%H(-2x2  +  Ww) 


(3-6) 


The  output  can  now  be  written  in  the  following  state  space  form: 


e 

'  i  M 

Xl 

t 

"  0  0 

'  5e  ' 

e 

0  -2/8J 

_  x2  . 

"T 

i  ./O' 

Ww 

The  second  augmentation  option  is  to  inject  the  driving  noise  at  the  aircraft  input 
as  shown  in  Figure  3-4. 


Figure  3-4.  Noise  Injected  at  Aircraft  Input 


Because  the  two  transfer  functions  act  in  series,  they  require  some  manipulation.  The 
state-space  equation  of  the  filter  remains: 


x2  =  -2x2  +  Ww  and  y-Wc-  J$Hx2 


(3-8) 


The  aircraft  input,  however,  is  now  5+  Wc,  and  Yc  becomes: 


xi  =  S«  +  Wc  and  y  =  e  =  x i 


(3-9) 


Combining  Equations  3-8  and  3-9  gives  the  state-space  representation  shown  below: 


*i 

*2 


0  JsJ 

Xl 

+ 

’  1  0  ' 

’  be  ' 

0  -2 

_  x2  . 

0  1 

Ww 

(3-10) 


To  find  the  output  vector  note  that 

e  =  Q  =  xi  and  e  =  6  =  x\  =  Jslsxz  +  8*  (3-11) 


Using  Equation  3-11,  the  output  vector  can  be  written  in  the  following  form: 


y  = 


e 

e 


1  0 
0  JsJ 


Xl 

*2 


de 

Ww 


(3-12) 


The  final  augmentation  option  requires  splitting  the  aircraft  into  two  separate 
transfer  functions,  one  before  and  one  after  the  noise  as  shown  in  Figure  3-5. 


Figure  3-5.  Noise  Injected  in  the  Middle  of  the  Aircraft  Dynamics 


This  form  is  not  practical  to  implement  for  the  integrator  example,  and  will  not  be 
considered  further  in  this  section. 
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The  linear  quadratic  regulator  (LQR)  problem  must  now  be  solved  so  that  the 
muscular  system  is  modeled.  This  is  accomplished  by  augmenting  the  aircraft  with 
integrators  and  iterating  the  control-rate  weighting  in  the  performance  index1.  For  this 
analysis,  the  noise  will  be  injected  at  the  aircraft  input,  and  Equations  3-10  and  3-12  will 


be  used.  The  noises  (including  W  )  can  be  set  to  zero  for  the  LQR  solution,  and  the 


resulting  state-space  is: 
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This  state  space  equation  can  be  augmented  with  integrators  by  adding  a  third 
state,  x3  =  5e,and  driving  the  system  with  5* .  The  resulting  state-space  is: 
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Equation  3-14  is  now  used  to  find  the  LQR  solution.  Because  output  weighting  is 
desired,  the  user  must  select  a  state  deviation  weighting  matrix,  Q,  where 

Q=CTQyC  (3-15) 

Qy  =  diag{q,  0}  (3-16) 

1  Consult  the  discussion  in  Appendix  A  for  further  information  on  this  method. 
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q  is  always  a  scalar  value.  In  this  example  q  =  1  will  be  used.  The  weighting  on  error 
rate  is  automatically  set  to  zero  in  this  model. 

Now  the  control  weighting  is  iterated  until  the  gain  on  the  last  state  equals  the 
inverse  of  the  desired  neuromuscular  time  constant.  The  optimal  control  law  from  the 
LQR  solution  of  Equation  3-14  is: 


U-be~  ~KCX 

(3-17) 

Partitioning  the  gain  and  state  matrices  results  in 

8«  =  -[  Kci  KC2  ]  * 

(3-18) 

8« =  ~Kci  •  x  —  Kc 2  •  be 

(3-19) 

K-el‘be+be  =  -K£Kcl  X 

(3-20) 

Equivalently,  the  following  relationship  exists  between  the  control  input  before  and  after 

the  neuromuscular  system: 

5e  _  1 

b'  e  X„S+\ 

(3-21) 

be  is  the  control  input  before  the  muscular  lag  and  5e  is  the  control  input  after  the  lag. 

Equation  3-21  can  be  rewritten  as: 

T*8e+8e  =5£ 

(3-22) 

The  terms  in  Equations  3-20  and  3-22  can  be  equated.  This  results  in  the  following 
relationships: 


t H=K£  and  8'  e  =  -K~c\  KciX  (3-23) 
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In  Equation  3-18,  Ka  is  the  gain  on  the  last  state.  Thus,  the  neuromuscular  lag  can  be 
modeled  by  iterating  the  weighting  on  control  rate  until  K  '^  =  Tn  (2:10). 

Continuing  with  the  integrator  example,  with  q  =  1  and  Tn  =  0.08,  the  STI  model 
will  return  a  value  for  g  of  0.00017.  Further  insight  can  be  gained  by  solving  this 
problem  using  the  MATLAB™  command,  Iqr.  Using  Equation  3-14  and  the  weighting: 


’  1  0  " 
0  0 

and  Q=CTQyC  = 

’  1  0  0  " 
0  0  0 

0  0  0 

yields  the  gain  matrix: 


K=[  76.70  31.03  12.39  ] 


(3-24) 


(3-25) 


The  inverse  of  the  last  gain  is  0.08,  the  desired  muscular  lag1.  The  corresponding 
closed-loop  transfer  functions  are: 


r 

r  (*+2)  1 

u 

(j+2X-J2+12.4j+76.7) 

£ 

j(j+2) 

L  u  J 

_  (j+2Xj2+12.4j+76.7) 

(3-26) 


The  complex  pole  has  a  natural  frequency  of  8.76  radians  per  second  and  represents  the 
neuromuscular  dynamics.  The  iterative  method  used  by  this  model  will  always  return  a 
pair  of  complex  poles  slightly  below  the  desired  neuromuscular  break  frequency. 

Compare  this  result  with  the  results  achieved  by  appending  the  neuromuscular 
dynamics  directly  to  the  aircraft  dynamics.  The  first  order  muscular  lag  can  be  written  as 
shown  in  Equation  3-27. 


1  The  STI  literature  refers  to  the  Time  Constant  (tn)  using  the  form:  (rNs  +1 ). 
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=>  X3=-—x+j^Si  and  6*  =  x3 


(3-27) 


5,  1 

8/e  Xn  S  +  1 


Placing  Equations  3-13  and  3-27  in  series  and  using  the  identity 


Pl=p2-Px  = 

A2  BiC\ 

0  Ax 

B2D\ 

Bx 

C2  D2Cx 

D2Dx 

(3-28) 


where 


‘ A , 

Bi  ' 

L  Ci 

Di . 

results  in  the  following  state-space  representation: 
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(3-29) 


Using  Equation  3-29,  Tn  =  0.08,  and  q  =  g  =  1,  the  LQR  solution  is: 


K=[  1  1.0563  .0770  ] 


(3-30) 


12.5(5+2) 

U 

(5+2)(j+1)(5+12.46) 

e 

12.5-5(5+2) 

L  u  J 

_  (5+2)(5+l)(5+12.46) 

(3-31) 


While  this  approach  places  a  pole  very  near  the  desired  muscular  lag  frequency  (12.5),  it 
allows  the  other  pole  to  float  well  inside  the  pilot  bandwidth.  The  iterative  method 
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produced  a  40  dB  per  decade  roll-off  slightly  below  the  desired  muscular  frequency. 

This  approach  produced  a  20  dB  per  decade  roll-off  at  the  desired  muscular  frequency. 
The  first  choice  best  models  the  human  pilot  and  is  used  in  the  STI  optimal  pilot  model. 

The  Kalman  filter  problem  must  now  be  solved.  Since  this  is  a  single  axis 
example,  the  fractional  attention  parameter,  f  is  one.  For  simplicity,  this  example  will 
assume  the  indifference  thresholds  are  zero.  The  noise  intensity  equations  from  Figure 
3-2  take  the  following  form: 

Vyi  =  diag{ Vy  ivy2}  (3-32) 

Vyi  =  pyi  •  n  •  <s2yi  and  Vu=puTtal  (3-33) 

In  these  equations  the  noise  ratios  (p)  are  user-defined,  constant  values.  These  values  are 
entered  in  dB,  but  it  should  be  noted  that  this  model  uses  power  spectra  dB  according  to 
the  following  relationship. 

dB=  10  •  log10(x)  (3-34) 

This  example  uses  the  experimentally  estimated  noise  ratios  of  0.01  and  0.003  (6:18),  or: 

Pyi  =  -20  dB  and  p„=-25  dB  (3-35) 

When  the  program  is  initiated,  an  initial  guess  for  the  noise  intensities  must  be  made. 

The  LQG  problem  is  solved  and  the  variances  of  the  states  are  computed.  The  noise 
intensities  are  then  iterated  until  Equations  3-33  is  satisfied  within  some  tolerance. 
Assuming  a  pilot  delay  of  0.15  seconds  and  unit-intensity  driving  noise  (Vw),  the 
optimal  pilot  model  converged  after  seven  iterations.  The  information  in  Figure  3-6  was 
displayed. 
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Iteration  #  7 

Optimal  cost:  q,  g,  output,  input  rate,  total 

1.0000E+00  1.7000E-04  1.1803E-01  4.1207E-02  1.5923E-01 

Performance:  E{y_lA2},  E{y_2A2},  E{u_aA2},  E{uA2},  E{(du/dt)A2} 

1.1803E-01  3.0834E+00  4.8469E+00  3.8633E+00  2.4240E+02 

Old  noise  intensities  (V_yl,V_y2,V_ua):  3.7167E-03  9.6967E-02  4.8176E-02 

New  noise  intensities  (V_yl,  V_y2,  V_ua):  3 .7079E-03  9.6869E-02  4.81 52E-02 

Noise  ratios  dB(Ao_yl,rho_y2,rho_ua):  -19.9897  -19.9956  -24.9978 

Max.  noise  ratio  difference  ■*  1. 02806 IE-02  dB,  Threshold  =  .1  dB 
Finished  with  iterations 

Figure  3-6.  STI  Optimal  Pilot  Model  Output  Display 

The  performance  index  value  in  the  top  right-hand  portion  of  the  display 
(1.5923E-01)  was  found  from  the  following  relationship: 

J=E{y\)-q  +  E{u*}g  (3-36) 

Note  that  by  converting  the  noise  ratios  (p)  back  into  decimal  form  (0.01, 0.01,  and 
0.001),  the  noise  intensities  can  be  verified  using  the  following  equations. 

Vyi  =£{yi}pi7t  Vy2  =E{y\}p2%  Vua  =  (3-37) 

Unfortunately,  due  to  the  way  the  display  code  was  written,  the  above  equations  only 
hold  when  the  Old  noise  intensities  line  is  used.  The  final  solution  and  performance 
index  value  are  unaffected  by  this  oddity,  however. 

This  model  requires  a  non-zero  pilot  delay.  It  finds  the  LQG  solution  in  the 
digital  domain  using  a  sampling  period  equal  to  the  pilot  delay.  If  the  delay  is  zero,  the 
numerical  algorithm  called  by  the  model  will  give  an  error  message  and  the  optimal  pilot 
model  will  display  erroneous  data. 
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The  STI  optimal  pilot  model  uses  the  compensator  gains  to  form  a 
transfer-function  representation  of  the  pilot.  This  is  accomplished  by  modeling  the  pilot 
delay  with  a  second  order  Pade  approximation. 

The  performance  index  can  be  used  to  estimate  an  Cooper-Harper  rating  using 
Equation  3-2.  In  this  equation,  the  rating  is  scaled  by  the  forcing  function  bandwidth, 
<dw,  as  well  as  the  forcing  function  root  mean  square  (RMS)  error  amplitude,  oc.  The 
RMS  value  is  found  by  noting  that  for  linear  systems  driven  by  white  noise  w(t): 


where 


x  =  Ax+Ew(t) 
y=Cx 


(3-38) 


«)}  =  o 

E{Mit)wT(t  +  x))  =  Q-8(x) 

w(t)  =  Zero  Mean  Stationary  White  Gaussian 
E{w(t)}=  Gaussian  Probability  Density 

Q  =  Constant  Intensity  Matrix  for  w(t) 


The  covariance  matrix  X  is  defined  as: 


X=E{x(t)xT(i)} 


(3-39) 


and  is  the  solution  of  the  following  Lyapunov  equation: 


A  X  +X -AT+  E  Q  E=0 


(3-40) 


The  output  covariance  matrix  is  then 


Y=E{y(t)yT(t)}=CX.CT 


(3-41) 


Because  the  STI  optimal  pilot  model  requires  single-input-single-output  dynamics,  Fis  a 
scalar  and  is  the  output  variance,  o2 ,  of  the  forcing  function  filter  ( 1 1 :97-l  1 1).  The 
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Program  CC  command  -  MEAN,  Yw  -  will  compute  this  value  as  long  as  the  task  forcing 
function  coloring  filter  is  driven  by  unit  intensity  white  noise  (Vw=  1). 

In  this  example  the  RMS,  oc,  is  1.48  and  the  filter  bandwidth,  ©w,  is  2.2  radians 
per  second.  Substituting  these  values  into  Equation  3-2  result  in  the  following  equation. 

Rating  =  1.7  +  3.7  •  log  l0(J)  (3-42) 

For  a  performance  index  value  of  0.159  (see  Figure  3-6),  the  predicted  Cooper-Harper 
rating  is  -1.25.  Obviously,  this  is  not  a  practical  handling  qualities  rating.  The  aircraft 
dynamics  and  tracking  task  forcing  function  used  for  this  example  were  overly 
simplified.  More  importantly,  the  parameters  used  when  running  this  model  were  chosen 
to  parallel  an  example  in  Reference  25.  As  shown  in  the  next  section,  these  parameters 
must  be  more  carefully  selected. 

Sensitivity  Analysis.  The  pilot  describing  function  and  Cooper-Harper  rating 
predicted  by  the  STI  optimal  pilot  model  are  sensitive  to  the  many  parameters  the  user 
must  select.  This  section  conducts  a  sensitivity  analysis  of  the  pilot  model  parameters 
and  offers  some  selection  guidance  in  an  effort  to  make  this  model  more  practicable. 

The  following  pitch-attitude  transfer  functions  were  used  to  conduct  the 
sensitivity  analysis. 


Case  1: 

0  100 
be  s(s+100) 

(3-43) 

Case  2: 

0  20(s  + 1 .25)  £_<m, 

8e  s(s2  +  8s  +  25) 

(3-44) 

Case  3: 

0  20(s+1.25) 

8e  s(s2+8s  +  25) 

(3-45) 
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Case  4:  9  20(5+1.25)  g.033f 

5«  s(s2  + 1.85  +  25) 

Case*  8,  20(5+ !. 25) 

,y(s2  +  1.8.J  +  25) 


(3-46) 

(3-47) 


These  equations  of  motion  were  chosen  for  two  reasons.  First,  they  represent  a  broad 
range  of  aircraft  dynamics  and  are  therefore  well  suited  for  the  sensitivity  analysis. 
Second,  these  dynamics  were  evaluated  in  a  previous  experiment  using  the  Large 
Amplitude  Multi-Mode  Aerospace  Research  Simulator  (LAMARS)  (Reference  19). 

All  of  the  dynamics  have  fighter  type  responses.  Case  1  is  nearly  an  integrator. 
The  other  cases  have  the  characteristics  shown  in  Table  3-1 . 


Table  3-1 

STI  Optimal  Pilot  Model  Evaluation  Dynamics 


Case  2 

Short  Period 
Damping  Ratio, 

Short  Period  Natural 
Frequency,  a>„ 

Delay,  x 

0.8 

5 

0.033 

Case  3 

0.8 

5 

0.2 

Case  4 

0.18 

5 

0.033 

Case  5 

0.18 

5 

0.2 

These  dynamics  were  evaluated  by  McRuer  using  a  sum-of-sines  tracking  task  generated 
on  the  LAMARS  head-up  display  (HUD),  and  were  assigned  the  Cooper-Harper  ratings 
shown  in  Table  3-2  (19:19). 


Table  3-2 

Cooper-Harper  Ratings  Used  for  Sensitivity  Analysis 


■BBT 

Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

jRating 

2 

3 

4 

4 

6 
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The  user  must  choose  the  following  parameters  before  running  the  STI  optimal 


pilot  model. 

1 .  Aircraft  Model  Order  ( Yc ) 

2.  Task  Forcing  Function  (Yw) 

3.  Intensity  of  the  Driving  Noise  (Vw) 

4.  Filter  Augmentation  Method  (Aircraft  Output  or  Input) 

5.  Neuromuscular  Time  Constant  (TJ 

6.  State  Deviation  Weighting  (q) 

7.  Initial  Guesses  and  Convergence  Thresholds  for  Control  Rate 
Weighting  and  Noise  Intensities  (g,  Vy ,  and  Vu) 

8.  Pilot  Delay  (x) 

9.  Observation  and  Motor  Noise  Ratios  (py, ,  py2 ,  and  pu ) 

10.  Visual  Indifference  Thresholds  (Tyl  and  Ty2) 

1 1 .  Fractional  Attention  Parameter  (f) 

The  data  gathered  in  this  thesis  indicates  that  the  values  in  Table  3-3  should  be  used. 


Table  3-3 

Recommended  Parameters  for  STI  Optimal  Pilot  Model 


Parameter 

Symbol 

Recommended  Value 

Aircraft  Model  Order 

Yc 

Lowest  Feasible 

Task  Forcing  Function 

Yw 

J2 

6.25 -s2  +3.54-5+1 

Intensity  of  the  Driving  Noise 

vw 

1 

Filter  Augmentation  Method 

Output 

Neuromuscular  Time  Constant 

0.08 

State  Deviation  Weighting 

q 

1 

Initial  Guesses  for  Control  Rate 
Weighting  and  Noises 

&.VV* 

ICTi 

fBAnssuHs 

Convergence  Thresholds  for  Control 
Rate  Weighting  and  Noises 

— 

g:  0.001 

Vv,Vu:  0.1 

Pilot  Delay 

X 

0.2  seconds 

Observation  and  Motor  Noise  Ratios 

-20  dB  for  All 

Visual  Indifference  Thresholds 

0 

Fractional  Attention  Parameter 

f 

1 

In  the  following  discussion  each  of  these  parameters  is  evaluated  individually  with  the 

others  held  constant  at  the  values  shown  in  this  table1. 

Aircraft  Model  Order.  The  lowest  model  order  consistent  with  the  task 

should  be  used.  High-order  dynamics  should  be  modeled  by  an  equivalent  delay  if 

possible  using  a  lower  order  equivalent  system2.  The  STI  optimal  pilot  model  computes 

a  predicted  pilot  describing  function  of  order  2n+5  where  n  is  the  number  of  aircraft  and 

filter  states.  For  example,  the  model  returned  a  fifteenth  order  pilot  describing  function 

when  the  short  period  approximation  of  Case  2,  was  analyzed. 

Task  Forcing  Function.  The  task  forcing  function  recommended  by  this 

thesis  is  a  second  order  Butterworth  filter,  shown  in  Equation  3-48,  with  a  break 
frequency  of  0.4  radians  per  second.  The  numerator  gain  (yj )  makes  the  filter 

bandwidth,  coM ,  equal  to  the  break  frequency. 

Yw  = - &= -  (3-48) 

.  j  . 2  J2  , 

U)  o>„  1 

Table  3-4  shows  the  relationship  between  filter  bandwidth  and  predicted 
Cooper-Harper  rating.  The  forcing  function  root  mean  square  (RMS)  errors  were 
computed  using  the  Program  CC  command  -  Mean,yw  -  and  the  predicted 
Cooper-Harper  ratings  were  computed  using  Equation  3-2.  The  pure  delays  in  Equations 
3-44  through  3-47  were  added  to  the  pilot  delay,  t.  When  modeled  instead  by  a  first 

order  Pade  approximation,  the  results  were  within  one-hundredth  of  a  rating. 

1  A  macro  containing  typical  Program  CC  commands  used  for  this  analysis  is  presented  at  the  end  of 
Appendix  B. 

2  The  lower  order  equivalent  system  match  routine  described  in  MIL-STD-1797  was  implemented  on  PC 
MatLab  in  a  Handling  Qualities  Toolbox  for  this  thesis.  Copies  of  this  toolbox  are  available  from  the 
author. 
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Table  3-4 

The  Effects  of  Forcing  Function  Bandwidth  on  Predicted  Cooper-Harper  Rating 

—  STI  Optimal  Pilot  Model 


Bandwidth, 

ow 

RMS 
Error,  ac 

Case  1 
(2)' 

Case  2 
(3)1 

Case  3 
(4)1 

Case  4 
(4)1 

Case  5 
(6)' 

Legend* 

0.1 

0.266 

3.99E-5 

1.17E-4 

1.86E-4 

3.50E-4 

4.33E-4 

J 

0.9 

2.6 

4.4 

WESM 

Rating 

0.4 

0.532 

4.73E-3 

9.72E-3 

1.66E-2 

2.33E-2 

3.09E-2 

J 

1.9 

3.0 

3.9 

4.4 

4 

.9 

Rating 

0.5 

0.595 

1.00E-2 

1.92E-2 

3.24E-2 

4.30E-2 

5.7E-2 

J 

2.0 

3.0 

3.9 

4.3 

4.8 

Rating 

0.8 

0.752 

4.62E-2 

7.45E-2 

1.23E-1 

1.41E-1 

1.86E-1 

J 

2.2 

3.0 

3.8 

4.0 

4.4 

Rating 

1 

0.841 

9.23E-2 

1.38E-1 

2.21E-1 

2.35E-1 

3.09E-1 

J 

2.2 

2.9 

3.6 

4 

.2 

Rating 

2 

1.19 

6.63E-1 

7.73E-1 

1.07 

9.35E-1 

1.15  ! 

J 

2.1 

2.3 

2.8 

2.6 

2.9 

Rating 

5 

1.88 

4.16 

3.76 

3.69 

3.61 

3.60 

J 

0.6 

0.4 

0.4 

0.4 

0.4 

Rating 

1  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 

2  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 


The  data  from  Table  3-4  is  plotted  in  Figure  3-7  for  clarity.  Close  examination  of 
this  figure  reveals  that  the  STI  optimal  pilot  model  responded  to  the  forcing  function 
characteristics  in  much  the  same  way  as  a  human  pilot  would.  As  the  bandwidth  of  the 
forcing  function  approached  the  pilot  bandwidth,  the  difference  between  the  predicted 
ratings  became  negligible.  Likewise,  the  human  pilot  will  not  attempt  to  follow  the 
commands  as  this  frequency  is  approached,  and  is  therefore  not  able  to  distinguish 
between  the  different  dynamics. 
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Figure  3-7.  The  Effects  of  Forcing  Function  Bandwidth  on  Predicted  Cooper-Harper 

Ratings  —  STI  Optimal  Pilot  Model 

The  logical  first  choice  for  forcing  function  bandwidth  is  one  that  mirrors  that  of 
the  actual  tracking  task.  The  STI  optimal  pilot  model  does  not  work  well  when  this 
choice  is  made.  McRuer's  evaluation  of  these  dynamics  used  a  forcing-function 
bandwidth  of  2  radians  per  second  which  is  typical  for  compensatory  tracking  tasks.  As 
shown  in  Figure  3-7  the  ratings  in  this  region  are  tightly  grouped  and  are  not  good 
predictors  of  the  experimental  ratings.  Instead,  the  best  results  will  be  obtained  if  each 
case  is  evaluated  for  a  range  of  forcing  function  bandwidths,  and  the  one  producing  the 
worst  Cooper-Harper  rating  is  used.  Note  that  in  Figure  3-7  this  maximum  occurs  over  a 
broad  frequency  making  the  top  of  each  curve  fairly  flat.  A  break  frequency  of  0.4 
radians  per  second  worked  well  for  all  five  cases  examined  in  this  chapter. 

Driving  Noise  Intensity.  The  intensity  of  the  driving  noise,  Vw,  can  be 
included  either  as  a  numerator  gain  in  the  task  forcing  function,  Yw ,  or  directly  by 
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varying  the  noise  intensity  value,  Vw.  Both  methods  yield  identical  results.  In  either 
case  the  predicted  rating  is  independent  of  the  noise  intensity. 

The  method  for  finding  the  forcing  function  root  mean  square  (RMS)  error 
amplitude,  oc,  described  in  Equations  3-38  through  3-41  and  the  Program  CC  command 
that  implements  this  method  —  Mean,  Yw  —  assume  a  unit  intensity  driving  noise,  Vw .  If 
the  driving  noise  is  not  unit  intensity,  the  numerator  of  the  task  forcing  function,  Yw, 
must  be  multiplied  by  the  square  root  of  the  driving  noise  intensity,  Vw ,  to  form  an 
equivalent  system.  It  is  therefore  easier  to  always  use  unit  intensity  driving  noise. 

Because  the  Cooper-Harper  rating  formula.  Equation  3-2,  is  normalized  for  ac , 
the  predicted  Cooper-Harper  rating  is  independent  of  the  numerator  gain  of  the  task 
forcing  function,  Yw .  Table  3-5  shows  the  predicted  Cooper-Harper  ratings  for  a  range  of 
numerator  values.  Note  that  while  the  performance  index,  J,  increased  with  the  forcing 
function  RMS  error  amplitude,  ac,  the  predicted  Cooper-Harper  ratings  did  not  change. 


Table  3-5 

The  Effects  of  Task  Forcing  Function  Numerator  on  Predicted  Cooper-Harper  Rating 

--  STI  Optimal  Pilot  Model 


Yw 

Numerator 

Case  1 
(2)2 

Case  2 
(3)2 

Case  3 
(4)2 

Case  4 
(4)2 

Case  5 
(6)2 

Legend* 

0.532 

4.73E-3 

9.72E-3 

1.66E-2 

2.33E-2 

3.09E-2 

J 

1.9 

mm 

3.9 

4.4 

4.9 

Rating 

2 

0.752 

9.43E-3 

1.94E-2 

3.32E-2 

4.68E-2 

6.17E-2 

J 

1.9 

3.0 

3.9 

4.4 

4.9 

Rating 

4 

1.50 

3.78E-2 

7.79E-2 

1.33E-1 

1.87E-1 

2.47E-1 

J 

1.9 

3.0 

3.9 

4.4 

4.9 

Rating 

'  Forcing  function  bandwidth  was  0.4  radians  per  second  for  all  runs  in  this  table. 

3  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 
3  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 
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Since  the  numerator  gain  of  the  task  forcing  function,  Yw ,  has  no  effect  on  the  predicted 
Cooper-Harper  rating,  the  square  root  of  two  was  used.  This  value  makes  the  bandwidth 
—  defined  as  the  zero  magnitude  (dB)  frequency  --  equal  to  the  filter  break  frequency. 

Filter  Augmentation  Method.  As  discussed  earlier  in  this  chapter,  the 
driving  noise  can  be  added  at  the  aircraft  input,  output,  or  between  to  parts  of  the  aircraft 
dynamics.  This  study  chose  output  augmentation  for  two  reasons.  First,  it  makes  better 
physical  sense.  Second,  the  Cooper-Harper  ratings  predicted  using  the  output 
augmentation  method  better  matched  the  experimental  results. 

The  STI  optimal  pilot  model  uses  only  error  and  error  rate  for  feedback.  When 
the  output  augmentation  method  of  Figure  3-3  (page  3-9)  is  chosen,  this  error  is  the 
difference  between  the  colored,  forcing  function  driving  noise  and  the  aircraft  output. 
This  is  similar  to  an  actual  tracking  task,  including  the  LAMARS  task,  where  the  pilot 
perceives  a  difference  between  his  actual  condition  and  his  desired  position  and  acts  to 
minimize  this  error.  If  the  noise  is  injected  at  the  aircraft  input  as  shown  in  Figure  3-4 
(page  3-10),  this  physical  significance  is  lost.  Injecting  the  noise  between  two  parts  of 
the  aircraft  dynamics,  as  shown  in  Figure  3-5  (page  3-11),  is  for  specialized  analysis  and 
not  considered  further  in  this  thesis. 

The  Cooper-Harper  ratings  predicted  using  the  output  augmentation  method 
better  matched  experimental  results.  The  predicted  Cooper-Harper  ratings  for  both  input 
and  output  noise  injection  are  shown  in  Table  3-6.  Notice  that  input  injection  values  are 
overly  pessimistic  when  compared  to  experimentally  obtained  values. 


3-26 


Table  3-6 

The  Effects  of  Noise  Injection  Method  on  Predicted  Cooper-Harper  Rating 

—  STI  Optimal  Pilot  Model 


Noise  Injection 
Method 

wBm 

mm 

B99I 

HI 

mm 

Legend1 

Input 

1.20E-2 

5.13E-2 

8.28E-2 

1.67E-1 

2.11E-1 

J 

3.3 

mm 

6.5 

7.6 

8.0 

Rating 

Output 

4.73E-3 

9.72E-3 

1.66E-2 

2.33E-2 

3.09E-2 

J 

1.9 

3.0 

3.9 

warn 

4.9 

Rating 

1  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 
1  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 


Neuromuscular  Time  Constant.  The  neuromuscular  time  constant  is  used 
to  limit  the  pilot  model  bandwidth.  The  actual  bandwidth  for  this  model  was  lower  than 
the  time  constant  choice  should  have  produced.  Thus,  a  time  constant,  Tn ,  of  0.08,  or  a 
muscular  lag  at  12.5  radians  per  second,  is  recommended. 

The  iterative  method  of  using  control  rate  weighting  to  model  the  neuromuscular 
lag  always  produces  a  pair  of  complex  poles  at  a  frequency  less  than  the  desired  muscular 
lag  frequency.  As  shown  with  the  integrator  example  in  Equation  3-26  (page  3-14),  a 
neuromuscular  time  constant,  Tn,  of  0.08  placed  a  pair  of  complex  poles  at  8.76  radians 
per  second.  The  actual  human  neuromuscular  time  constant  has  been  experimentally 
measured  between  0.08  and  0.12,  resulting  in  a  first  order  lag  between  8.3  and  12.5 
radians  per  second  (16:29).  The  high  end  of  this  range  (Tn=  0.08)  was  chosen. 

State  Deviation  Weighting.  It  makes  no  difference  what  value  is  used  for 
the  state  deviation  weighting.  The  STI  model  results  are  independent  of  this  value. 

Because  the  STI  optimal  pilot  model  requires  single-input-single-output  aircraft 
dynamics,  the  state  deviation  weighting  matrix  is  formed  by  selecting  a  scalar  value,  q,  as 
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shown  in  Equation  3-16  (page  3-12).  The  model  then  finds  the  control  rate  weighting,  g, 
that  establishes  the  desired  neuromuscular  lag.  The  ratio  between  the  control-rate  and 
state  deviation  weighting  is  fixed  by  this  process. 

The  gain  matrix  shown  in  Equation  3-25  (page  3-14),  was  computed  in  the 
integrator  example  presented  earlier  in  this  chapter.  In  this  example,  the  state  deviation 
weighting,  q,  was  one  and  the  state  deviation  weighting,  g,  used  to  model  the 
neuromuscular  lag  was  0.00017.  When  the  same  problem  was  repeated  with  a  state 
deviation  weighting  of  5,  the  model  set  the  control  rate  penalty,  g,  at  0.00085  to  model 
the  muscular  dynamics  and  computed  the  same  gain  matrix.  This  weighting  was  five 
times  the  previous  value,  just  as  the  state  deviation  penalty  was  five  times  that  in  the 
previous  example. 

Initial  Guesses  and  Convergence  Thresholds.  Because  the  STI  optimal 
pilot  model  finds  the  control  rate  weighting,  g,  and  the  noise,  and  Vu ,  iteratively,  the 
model  requires  initial  guesses  and  convergence  thresholds  for  each.  Within  reason,  it 
does  not  matter  what  values  are  chosen  for  the  initial  guesses.  The  number  of  iterations 
required  for  convergence  does  not  increase  significantly  if  a  poor  initial  guess  is  made. 
This  thesis  used  the  convergence  thresholds  recommended  in  Reference  25  (25:29). 

Pilot  Delay.  The  delay  time  of  the  human  pilot  has  been  measured 
experimentally  between  0. 1 5  and  0.3  seconds  (6: 1 8).  A  pilot  delay  of  0.2  seconds  was 
used  in  this  analysis. 

The  effects  of  pilot  delay  on  predicted  Cooper-Harper  rating  are  shown  in 


Table  3-7. 
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Table  3-7 

The  Effects  of  Pilot  Delay  on  Predicted  Cooper-Harper  Rating 
—  STI  Optimal  Pilot  Model 


Pilot  Delay 
(seconds) 

0.10 

0.15 

0.20 

0.25 

0.30 

Legend2 

Case  1 
(2)' 

2.53E-3 

3.52E-3 

4.73E-3 

6.15E-3 

7.80E-3 

J 

0.9 

1.4 

1.9 

2.3 

MSM 

Rating 

Case  2 
(3)' 

6.28E-3 

7.93E-3 

9.72E-3 

1.16E-2 

1.37E-2 

J 

2.3 

3.0 

3.3 

3.6 

Rating 

Case  4 
(4)' 

1.75E-2 

2.06E-2 

2.33E-2 

2.59E-2 

2.82E-2 

J 

4.0 

4.2 

4.4 

4.6 

4.7 

Rating 

1  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 
1  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 


Cases  3  and  5  are  not  included  in  this  table  because  their  dynamics  are  identical  to  those 
of  Cases  2  and  4.  Note  that  the  relationship  is  fairly  linear  with  every  ten  milliseconds  of 
added  delay  producing  no  more  than  a  tenth  of  a  rating  increase.  Also  note  that  added 
delay  does  not  affect  the  ratings  of  the  poor  dynamics  as  much  as  good  dynamics.  For 
the  same  added  delay,  the  predicted  rating  for  Case  1  dropped  from  0.9  to  2.7  while  the 
predicted  rating  for  Case  4  went  from  4.0  to  4.7. 

Noise  Ratios.  The  noise  ratios  are  perhaps  the  most  difficult  pilot  model 
parameters  to  select.  The  predicted  Cooper-  -atings  and  estimated  pilot  describing 
functions  are  especially  sensitive  to  these  values.  For  this  reason  it  is  desirable  to  stay  as 
close  as  possible  to  experimentally  measured  values.  Reference  6  recommended  using 
the  following  values  (6: 18). 

Pyi  =  -20  dB  and  p„  =  -25  dB  (3-49) 


The  STI  optimal  pilot  model  results  in  Table  3-8  were  found  using  these  values. 
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Table  3-8 

The  Effects  of  Noise  Ratios  on  Predicted  Cooper-Harper  Rating 
—  STI  Optimal  Pilot  Model  (pu  =  -25  dB) 


Bandwidth, 

<°w 

RMS 
Error,  oc 

Case  1 
(2)' 

Case  2 
(3)‘ 

Case  3 
(4)’ 

Case  4 
(4)' 

Case  5 
(6)' 

Legend1 

0.1 

0.266 

1.87E-5 

4.99E-5 

9.04E-5 

1.41E-4 

2.00E-4 

J 

— — 

1.2 

2.3 

2.9 

3.5 

Rating 

0.4 

0.532 

2.79E-3 

5.35E-3 

1.02E-2 

1.13E-2 

1.75E-2 

J 

1.0 

2.1 

3.1 

3.3 

4.0 

Rating 

0.5 

0.595 

6.11E-3 

1.11E-2 

2.11E-2 

2.11E-2 

3.38E-2 

J 

1.2 

2.2 

3.2 

3.2 

4.0 

Rating 

0.8 

0.752 

3.04E-2 

4.88E-2 

8.80E-2 

7.99E-2 

1.23E-1 

J 

1.5 

2.3 

3.2 

3.1 

3.8 

Rating 

1 

0.841 

6.31E-2 

9.55E-2 

1.67E-1 

1.42E-1 

2.16E-1 

J 

1.6 

2.0 

2.6 

2.2 

2.7 

Rating 

2 

1.19 

5.00E-1 

6.30E-1 

9.20E-1 

7.20E-1 

9.80E-1 

J 

1.6 

2.0 

2.6 

2.2 

2.7 

Rating 

5 

1.88 

3.50 

3.50 

3.60 

3.50 

3.50 

J 

0.3 

0.3 

0.4 

0.3 

0.3 

Rating 

1  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 

2  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 


As  in  Table  3-4,  the  results  are  displayed  for  range  of  forcing  function  bandwidths.  The 
data  from  Table  3-8  is  plotted  in  Figure  3-8  on  the  following  page  for  clarity.  In  all  cases 
the  predicted  Cooper-Harper  ratings  were  overly  optimistic  when  compared  with 
experimentally  obtained  ratings. 

Increasing  the  motor  noise  ratio,  pu ,  from  -25  dB  to  -20  dB  improved  the  model's 
predictions.  The  results  in  Table  3-4  and  Figure  3-7  (page  3-24)  were  obtained  with  a 
motor  noise  ratio  of  -20  dB,  and  were  more  in  line  with  the  experimentally  obtained 
ratings.  The  motor  noise  ratio,  pu,  was  used  in  this  thesis  to  fine  tune  the  STI  optimal 
pilot  model  predictions,  and  -20  dB  seemed  to  do  this  well. 


3-30 


6 


Figure  3-8.  The  Effects  of  Noise  Ratios  on  Predicted  Cooper-Harper  Ratings 
--  STI  Optimal  Pilot  Model  (pu  =  -25  dB) 

Visual  Indifference  Thresholds.  Visual  indifference  thresholds  in  the  STI 
optimal  pilot  model  tended  to  excessively  penalize  good  dynamics,  and  this  thesis 
recommends  using  thresholds  of  zero  unless  display  effects  are  being  specifically  studied. 

The  visual  indifference  thresholds  of  the  human  pilot  are  dependent  on  the  task 
display.  Dillow  and  Picha  used  0.05  degrees  and  0.18  degrees  per  second  as  typical 
values  for  pitch  axis  error  and  error  rate,  visual  indifference  thresholds  (6: 18).  The  STI 
optimal  pilot  model  results  found  using  these  values  are  presented  in  Table  3-9.  As 
shown  in  this  table,  the  predictions  obtained  using  indifference  thresholds  of  zero  better 
matched  experimental  results. 
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Table  3-9 

The  Effects  of  Visual  Indifference  Thresholds  on  Predicted  Cooper-Harper  Rating 

—  STI  Optimal  Pilot  Model 


Visual  Indifference 
Thresholds 

mm 

ISJIIJ 

w$m 

mSM 

mm 

wnm 

Legend1 

Error,  Tyl  =  0 

Error  Rate,  Ty2  =0 

4.73E-3 

9.72E-3 

1.66E-2 

2.33E-2 

3.09E-2 

J 

1.9 

3 

3.9 

4.4 

4.9 

Rating 

Error,  Tyl  =  0.05 
Error  Rate,  Ty2  =0.18 

1.14E-2 

1.60E-2 

2.32E-2 

2.78E-2 

3.51E-2 

J 

3.3 

3.8 

mm 

mm 

5.1 

Rating 

1  Cooper-Harper  rating  assigned  during  simulator  tracking  task  in  Reference  24 

2  Where  J  is  the  performance  index  and  Rating  is  the  predicted  Cooper-Harper  rating 


Fractional  Attention  Parameter.  The  fractional  attention  parameter  is  used 
to  predict  multi-axis  Cooper-Harper  ratings.  Each  axis  is  assigned  a  fraction  of  the  pilot's 
attention  so  that  the  fractions  total  one.  Because  the  STI  optimal  pilot  model  requires 
single  axis  dynamics,  the  fractional  attention  parameters  must  be  manually  iterated. 

In  a  two  axis  problem,  for  example,  the  user  must  find  the  performance  index  in 
each  axis  for  a  range  of  fractional  attention  parameters.  The  performance  indices  for 
those  combinations  of  fractional  attention  parameters  that  total  one  must  then  be 
summed.  The  lowest  total  performance  index  is  then  used  to  find  the  predicted 
multi-axis  rating  using  Equation  3-2.  Obviously,  this  procedure  becomes  more 
complicated  and  lengthy  as  additional  axes  are  added. 

Pilot  Describing  Function.  The  STI  optimal  pilot  model  computes  a  pilot 
describing  function  along  with  a  predicted  Cooper-Harper  rating.  This  transfer  function 
enhances  the  model's  abilities  as  a  predictive  tool.  Unfortunately,  these  transfer  functions 
are  very  high  order  and  often  do  not  reduce  to  the  classical  lead-lag-delay  pilot  model 
with  acceptable  accuracy. 
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Using  the  values  in  Table  3-3,  the  STI  optimal  pilot  model  predicted  the 
following  pilot  describing  functions  for  Cases  1  through  5  (YpI-  Yp5)\ 


_  174(0X.0448)[.707,.4](1.65X5.79X12.63)(12.91)[-.866, 17.3](100)2 
p  1  ~  (0)[.707,  .4]2(6.02X12.57X12.63)[.274, 18.961(45.4X94X100) 


Y„  2  = 


92.6(0)(,032)[.707,  .4](1.45)[.814,4.S61[.8, 5](5.08X12.5)2[-.866, 14.9] 
(0)[.707,  .4]2(1.25)(4.96)[.8, 5](12.43X12.46)[.146, 14.2J[.84,27.3] 


(3-50) 


(3-51) 


_92.7(0X.035)[.707,.41(U9)(4.1)[.79  6,4.781[.8,5][-.866,8.66](12.5)2 
"3  (0)[.707,  .4]2(1.25)(4.02)[.8, 5][.043, 1 1.5](12.36)(12.46)[.81, 23.5]  1 

_  46.2(0)(.028)[.707,.4](1.18)[.06,4.4](4.56X.i8,5](12.5)(12.6)[-.866,14.9] 
p 4  (0)[.707,  .4]2(1.25)(4.6)[.18,5](12.49X12.52)[.157, 12.8][.846, 23.3]  (  ’ 


Y  _  40,4(0X.028)[.707,.4](1.04X3.74)[.02S>4.96],[.18>5]>[-.866>8.66](12.S)2 
"5  (0)[.707,.4]2(1.25)(3.77)[.182,5](.082,10.7](12.4X12.5)[.8I8,19.12]  1 

The  pilot  delay  was  modeled  with  a  second  order  Pade  approximation.  There  are  several 

nea*-  pole-zero  cancellations  in  the  above  transfer  functions.  Thompson  also  recommends 

grouping  the  high-frequency  poles,  those  above  5  radians  per  second,  into  an  equivalent 

delay  (25:38).  An  alternative  approach  is  to  use  a  lower  order  equivalent  system  match 

routine  similar  to  that  used  in  MIL-STD-1797A  for  aircraft  dynamics  (Reference  5).  The 

following  equivalent  systems  were  computed  by  finding  the  best  fit  to  the  classical  gain, 


lead,  lag,  and  delay  pilot  model  between  0. 1  and  10  radians  per  second 2. 


3.756(s  + 2.921)  154j 

pl  (s+ 1.546) 


(3-55) 


1  Where  the  brackets  denote  [£,  toj  as  in  s1  +  2£ g>„+  a>J  ,  and  the  parentheses  denote  (t)  as  in  s  +  x. 

1  This  lower  order  equivalent  system  match  routine  was  implemented  on  PC  MatLab  in  a  Handling 
Qualities  Toolbox  for  this  thesis  and  is  available  from  the  author. 
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(3-56) 


1.29(5  +  3.248)  106j 

pl  (5  +  0.973) 


1.294(5  +  2.253) 
p3  (5  +  1.036) 


g-.l98s 


0.347(5  +  5.42)  OI23, 

p 4  (5  +  0.682) 


(3-57) 

(3-58) 


0.295(5  +  5.31)  13, 

pS  (5  +  0.765) 


(3-59) 


The  Bode  plots  of  the  predicted  pilot  describing  function  and  their  lower  order 
equivalent  systems  matches  are  shown  in  Appendix  B.  These  matches  were  especially 
poor  for  the  lightly  damped  dynamics  (Cases  4  and  5).  For  these  cases,  shown  in  Figures 
B-4  and  B-5,  the  STI  optimal  pilot  model  predicted  a  pilot  describing  function  that 
essentially  notched  out  the  aircraft  short  period  dynamics.  This  behavior  is  not  possible 
using  only  a  gain,  lead,  lag,  and  delay  model,  and  is  inconsistent  with  McRuer’s 
experimental  observations  (14:218)'. 


Summary 

The  optimal  pilot  model  was  developed  in  an  attempt  to  overcome  deficiencies  in 
the  classical  pilot  model  approach.  This  chapter  focused  on  an  optimal  pilot  model 
developed  by  Systems  Technology,  Incorporated  (STI),  for  use  with  Program  CC 
(Reference  25).  This  model  incorporates  nearly  every  important  aspect  of  other  optimal 
pilot  models  and  can  be  implemented  on  the  personal  computer,  but  it  lacks  parameter 
selection  guidance. 

1  Consult  the  discussion  on  page  2-13. 
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This  model  allows  the  feedback  of  only  one  variable  and  its  derivative,  therefore 
single-input-single-output  dynamics  must  be  used  for  both  the  aircraft  and  forcing 
function  dynamics.  The  model  solves  the  linear  quadratic  regulator  (LQR)  problem  by 
adjusting  the  control  rate  weighting,  g,  so  that  the  neuromuscular  dynamics  are  modeled. 
The  model  then  iteratively  solves  the  Kalman  filter  problem  so  that  the  desired  noise 
ratios  are  obtained.  Finally,  a  predicted  pilot  describing  function  is  computed  and  the 
performance  index  is  used  to  predict  a  Cooper-Harper  rating. 

The  user  must  choose  several  parameters  before  running  the  ST1  optimal  pilot 
model.  A  sensitivity  analysis  of  these  parameters  was  conducted  using  five  different 
pitch  axis  dynamics.  Based  on  this  analysis,  the  parameters  in  Table  3-3  (page  3-21) 
should  be  used  when  running  the  STI  optimal  pilot  model.  Additionally,  the  following 
conclusions  were  reached. 

1 .  Include  aircraft  delays  in  the  pilot  delay  or  model  them  with  a  first  order 
Pade  approximation.  The  results  are  the  same  for  both  methods  (page  3-22)1. 

2.  Model  the  task  forcing  function,  Yw ,  as  a  second  order  Butterworth  filter. 
Evaluate  the  dynamics  for  a  range  of  forcing  function  bandwidths,  and  use 
the  one  that  predicts  the  worst  Cooper-Harper  rating.  A  bandwidth  of  0.4 
radians  per  second  worked  well  for  all  five  cases  examined  in  this  chapter 
(page  3-24). 

3.  The  predicted  Cooper-Harper  ratings  are  independent  of  the  driving  noise 
intensity,  Vw  (page  3-25). 

4.  Injecting  the  forcing  function  driving  noise  at  the  aircraft  output  makes  better 
physical  sense  and  produces  more  accurate  predictions  (page  3-26). 

5.  Use  a  neuromuscular  time  constant,  Tn ,  of  0.08,  producing  a  lag  at  12.5 
radians  per  second.  The  optimal  pilot  model  establishes  the  model 
bandwidth  slightly  lower  than  it  should,  and  12.5  is  at  the  high  end  of  the 
experimentally  measured  range  (page  3-27). 


1  Refers  to  the  page  in  this  thesis  containing  this  conclusion. 
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6.  The  STI  optimal  pilot  model  results  are  independent  of  the  state  deviation 
weighting,  q  (page  3-28). 

7.  The  STI  optimal  pilot  model  is  insensitive  to  initial  guess  values  for  the 
control-rate  weighting  and  the  noise  intensities  (page  3-28). 

8.  Use  -20  dB  for  all  noise  intensities,  pyl ,  py2 ,  and  pu.  The  motor  noise 
intensity,  pu ,  was  increased  from  its  experimentally  measured  value  (-25  dB) 
to  fine  tune  the  model's  predictions  (page  3-30). 

9.  Set  the  visual  indifference  thresholds,  Tyl  and  Ty2 ,  to  zero  unless  display 
effects  are  being  studied.  Non-zero  values  overly  penalize  good  dynamics 
(page  3-31). 

10.  The  predicted  pilot  describing  functions  are  not  consistent  with  the  classical 
pilot  model  form  or  the  experimental  observations  of  McRuer  in 
Reference  14  (page  3-34). 
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4.  Flight  Test 


General 

As  shown  in  the  previous  chapter,  the  high  order  pilot  transfer  functions  predicted 
by  the  STI  optimal  pilot  model  were  not  consistent  with  classical  pilot  model  theory.  In 
an  effort  to  resolve  this  and  other  important  pilot  modeling  issues,  a  limited  evaluation  of 
human  pilot  response  was  sponsored  in  support  of  this  thesis  by  the  Air  Force  Flight 
Dynamics  Directorate.  This  flight  test  is  detailed  in  AFFTC-TLR-93-4I  (Reference  7). 
The  applicable  aspects  of  the  test  are  presented  in  this  chapter. 

This  chapter  is  divided  into  three  sections.  The  first  section.  Test  Procedures, 
provides  a  description  of  the  flight  test.  The  second  section.  Test  Results,  emphasizes 
those  results  applicable  to  optimal  pilot  modeling  theory.  This  section  focuses  on  single 
axis  tasks  and  contains  statistical  and  frequency  response  analysis  not  published  in 
AFFTC-TLR-93-41 .  Finally,  the  important  conclusions  of  this  chapter  are  summarized. 

Test  Procedures 

The  flight  test  was  conducted  at  Calspan  Corporation,  Buffalo,  New  York  by  a 
five  member  team  from  the  USAF  Test  Pilot  School  between  8  and  1 1  October  1993. 

Five  sorties  totaling  7.6  hours  were  flown  in  the  Calspan  variable  stability  Lear  II 
aircraft.  Ground  simulations  in  Lear  II  were  also  performed.  Four  different  pitch  and 
four  different  roll  axis  dynamics  were  evaluated  using  three  different  tracking  tasks.  For 
each  set  of  dynamics,  primary  pilot  response  parameters  were  recorded  and  examined 
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using  Fourier  transform  analysis  in  an  attempt  to  provide  an  insight  into  human  pilot 
behavior.  Pilot  comments  and  Cooper-Harper  ratings  were  also  recorded.  The  flight  test 
data  gathered  during  this  project  are  maintained  at  the  Air  Force  Flight  Dynamics 
Directorate  (WL/FIGC),  Wright-Patterson  AFB,  OH  45433,  and  are  available  for 
research  purposes.  AFFTC-TLR-93-41  serves  as  a  guide  for  the  flight  test  data  base  and 
provides  an  initial  look  at  the  test  results. 

Research  Vehicle  Description.  The  pilot  was  the  true  test  item  for  this  test 
program.  The  three  test  pilots  used  for  this  evaluation  had  a  variety  of  operational 
backgrounds.  Two  of  the  pilots  had  multi-engine  backgrounds.  One  had  extensive 
K.C-10  experience  and  the  other  had  extensive  C-130  experience.  The  third  pilot  was  an 
experienced  test  pilot  with  considerable  F-15  and  F-16  experience. 

The  Calspan  Variable-Stability  Lear  II  was  used  as  the  research  vehicle  for  this 
evaluation.  It  was  a  production  Leaijet  25  aircraft  that  was  extensively  modified  for  use 
as  an  in-flight  simulator.  The  basic  aircraft  is  shown  in  Figure  4-1 . 

The  safety  pilot's  (left  seat)  control  column  was  connected  directly  to  the  aircraft's 
control  surfaces  through  the  production,  reversible  push  rod  system  and  mirrored  the 
surface  positions.  The  evaluation  pilot's  (right  seat)  controls  were  removed  and  replaced 
with  a  variable,  center-stick  feel  system  that  was  part  of  the  in-flight  simulation  system. 

A  digital  computer  system  was  located  in  the  main  cabin.  It  was  designed  to 
model  aircraft  and  artificial  feel  systems  and  record  in-flight  data.  Real  time  monitoring 
of  up  to  64  selected  parameters  at  a  sampling  rate  of  100  hertz  was  possible.  The  data 
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Figure  4-1.  Two  Plan  View  of  Lear  25B 


were  stored  on  90  megabyte  personal  computer  (PC)  compatible,  removable  Bernoulli 
drive  cartridges  in  a  MATLAB™  compatible  format. 


A  color  flat  panel  display,  located  on  the  main  instrument  panel  in  front  of  the 
evaluation  pilot's  station,  displayed  the  tracking  tasks.  The  display  is  shown  in 
Figure  4-2. 


1 

0.5  deg.  Pitch  -*■  r 
5.0  deg.  Roll  * 


© 


Q  Command  Bar  for  Tacking  Tasks 
@  Fixed  Pipper  Reference 
®  Horizon  Line  (Extended  to  Width  of  Display) 


Figure  4-2.  Tracking  Task  Display 


Pitch  and  roll  errors  were  indicated  to  the  pilot  by  the  angular  deviation  between  the 
command  bar  and  the  extended  fixed-pipper  attitude  reference.  The  lengths  of  the 
subtends  on  the  attitude  reference  corresponded  to  0.5  degree  of  pitch  error  and  5  degrees 
of  roll  error.  A  horizon  line  was  also  displayed  to  reduce  the  potential  for  pilot 
disorientation.  A  0.025  second  delay  arose  from  the  display  process.  A  detailed 
description  of  the  Lear  II  and  the  flat  panel  display  is  given  in  AFFTC-TLR-93-41 
(7:33). 

Tracking  Tasks.  Three  different  types  of  tracking  tasks  were  used  for  this  test,  a 
discrete  tracking  task,  a  sum-of-sines  tracking  task,  and  a  regulator  task.  To  lessen  the 
pilot's  ability  to  memorize  the  tasks,  two  different  tasks  of  each  type  were  used  during 
testing.  Each  task  was  53  seconds  long. 
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Discrete  Tracking  Task.  The  discrete  tracking  task  consisted  of  a  series  of 
steps  and  ramps.  A  representative  task  is  shown  in  Figure  4-3.  The  initial  command  in 
each  axis  was  a  step.  The  maximum  commanded  input  was  ±2.25  degrees  from  the 
initial  condition  in  pitch  and  ±35  degrees  in  roll. 


Sum-of-Sines  Tracking  Task.  The  sum-of-sines  tracking  task  was  a 
random  appearing,  frequency-based  function  computed  using  the  following  formula: 

n 

Command  =  k  £a,  -  sin(© ,  +  4>,)  (4- 1 ) 

*=i 

The  tasks  used  in  this  test  were  formed  by  summing  13  sine  waves.  The  phases  (<{>;)  were 
randomly  chosen  and  the  task  gain  (K)  was  set  to  achieve  the  desired  task  amplitude. 

The  frequencies  (©;)  were  evenly  spaced  between  0.1  and  30  radians  per  second.  The 
amplitudes  (A;)  were  selected  using  a  comer  frequency  of  2  radians  per  second  and  a 
second-order  roll  off  producing  the  power  spectral  density  magnitudes  shown  in 
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Figure  4-4.  A  typical  task  is  shown  in  Figure  4-5.  While  this  task  appeared  random,  the 
power  spectral  density  was  concentrated  only  at  the  selected  frequencies  (cOj )  as  shown  in 
Figure  4-4: 


Figure  4-4.  Power  Spectral  Density  of  Sum-of-Sines  Tracking  Task 


Time,  t  (seconds) 


Figure  4-5.  Sum-of-Sines  Pitch  Axis  Tracking  Task 
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Regulator  Task.  The  regulator  task  was  computed  in  the  same  manner  as 
the  sum-of-sines  task.  Instead  of  driving  a  command  bar,  however,  this  task  was  input  as 
an  additive  to  the  pilot’s  stick  command.  It  had  the  effect  of  simulating  turbulence.  The 
pilot's  objective  during  this  task  was  to  maintain  wings  level,  zero  pitch  flight.  This 
tracking  task  had  the  same  frequency  content  as  the  sum-of-sines  task. 

Desired  And  Adequate  Criteria.  Each  airborne  test  case  was  examined  using  all 
three  types  of  tracking  tasks.  Only  the  discrete  and  sum-of-sines  tracking  tasks  were 
used  during  ground  simulation  runs.  Cooper-Harper  ratings  were  assigned  for  each  run 
in  accordance  with  the  following  criteria: 


Desired  Criteria:  Pilot  induced  oscillations  (PIO)  tendencies  do  not 
compromise  tracking  task.  Commanded  attitude  maintained  within  0.5 
degrees  in  pitch  and  5  degrees  in  bank  (measured  at  end  of  command  bar) 
for  50  percent  of  the  time  except  immediately  following  step  command 
changes.  See  Figure  4-6. 

Adequate  Criteria:  Commanded  attitude  maintained  within  1  degrees  in 
pitch  and  10  degrees  in  bank  (measured  at  end  of  command  bar)  for  50 
percent  of  the  time  except  immediately  following  step  command  change. 


Pitch 


Roll 


Within 
0.5  deg.0 


A 

T 


4 


Figure  4-6.  Tracking  Task  Desired  Criteria 


Dynamics.  The  following  dynamics  were  selected  based  on  past  Calspan  flight 
test  experience  so  that  Cases  1  and  A  would  produce  a  level  1  aircraft.  Subsequent  cases 
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degraded  either  the  damping,  roll  mode  time  constant,  or  system  delay  to  produce  the 
desired  spread  in  aircraft  handling  qualities. 


Longitudinal  Dynamics:  Lateral  Dynamics: 


Case  1: 

20(s+1.8)e-o-04j 
s(s2  +  8.45  +  36) 

Case  A: 

2.5e-°04 1 

5(5 +  2.5) 

(4-2) 

Case  2: 

20(s+  Lgy-0  04* 
s(s2  +  4.85  +  36) 

Case  B: 

e-0.04 s 

5(5+1) 

(4-3) 

Case  3: 

20(5  + 1 .8)e~°-24s 
s(s2  +  8.45  +  36) 

Case  C: 

2.5e~°24s 

5(5 +  2.5) 

(4-4) 

Case  4: 

20(5  + 1  .S)e~°-2is 
s(s2  +4.85  +  36) 

Case  D: 

e-0.24 j 

5(5+1) 

(4-5) 

These  cases  have  the  characteristics  summarized  in  Table  4-1. 

Table  4-1 

Case  Definition  Table 


Case  1 

Longitudinal 

1  Lateral 

Short  Period 
Damping,  C, 

Delay,  xe 

Roll  Mode  Time 
Constant,  TR 

Delay,  xe 

0.7 

0.04 

Case  2 

0.4 

0.04 

Case  3 

0.7 

0.24 

Case  4 

0.4 

0.24 

Case  A 

0.4 

0.04 

Case  B 

1.0 

0.04 

Case  C 

0.4 

0.24 

Case  D 

1.0 

0.24 
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These  dynamics  were  implemented  as  a  stick  position  command  system  as  described  in 
Appendix  C. 

During  the  simulations,  sideslip  angle  rate  (p)  was  driven  to  zero  by  the 
simulation  computers.  Additionally,  the  roll  and  yaw  axes  were  de-coupled  so  the  pilot 
could  fly  the  simulations  with  feet  on  the  floor. 

Data  Collection.  Both  ground  and  airborne  evaluations  were  conducted.  For  the 
ground  evaluations,  the  Lear  II  was  operated  in  the  ground  simulation  mode  inside  the 
Calspan  hanger.  In-flight  testing  was  conducted  under  day,  visual  meteorological 
conditions  (VMC)  with  no  more  than  occasional  light  turbulence. 

All  single  axis  ground  simulation  cases  were  flown  twice.  All  of  the  airborne 
single  axis  cases  were  evaluated  three  times.  In  planning  the  specific  test  points,  project 
engineers  ensured  adequate  pilot  variability  checks  by  assigning  some  pilots  the  same 
cases  twice  and  by  assigning  some  cases  to  more  than  one  pilot.  The  dynamics  cases 
were  also  evaluated  in  a  random  order. 

For  each  case,  the  test  pilots  flew  the  discrete  tracking  task,  followed  by  the 
sum-of-sines  tracking  task  and  the  regulator  task.  Each  task  was  53  seconds  long.  The 
regulator  tasks  were  only  evaluated  during  airborne  simulations  since  these  tasks  were 
designed  to  evaluate  normal  acceleration  feedback  cues. 

Immediately  after  each  task,  the  pilot  and  flight  test  engineer  completed  the  test 
point  comment  card  shown  in  Figure  C-3.  The  pilots  then  assigned  a  pilot  induced 
oscillation  (PIO)  rating,  using  the  scale  in  Figure  C-4,  and  a  Cooper-Harper  rating  using 
the  scale  in  Figure  2-1. 
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Twenty-two  parameters  were  recorded  with  a  sampling  rate  of  100  hertz  for  each 
test  case.  These  parameters  included  commanded  and  actual  pitch  and  roll  angle  as  well 
as  longitudinal  and  lateral  stick  force  and  deflection  (7:3-4).  The  percentage  of  time  the 
desired  and  adequate  criteria  were  met  as  computed  by  the  Lear  II  simulation  computer 
were  saved  in  separate  scoring  files. 

Data  Reduction .  In  an  effort  to  gain  insight  into  human  pilot  response,  a 
frequency  response  analysis  of  the  relationship  between  the  pilot’s  input  (task  error)  and 
the  pilot  output  (stick  force  and  displacement)  was  conducted.  This  analysis  was 
accomplished  using  the  MATLAB™  fast  Fourier  transform  software  described  in 
AFFTC-TLR-93-41  (7:55-58).  This  software  converted  time  domain  data  to  frequency 
domain  data.  In  doing  so,  it  yielded  the  power  spectral  densities  of  the  input  and  output 
parameters,  a  Bode  plot  of  the  transfer  function,  and  a  coherence  function.  This 
coherence  function  provided  a  measure  of  how  much  of  the  output  at  each  frequency  was 
caused  by  the  input  rather  than  by  noise  or  other  inputs.  A  coherence  value  of  1 .0  meant 
that  the  output  was  completely  a  function  of  the  input,  whereas  as  value  of  0.0  meant  that 
the  output  was  not  a  function  of  the  input.  This  software  was  included  in  the  data  base 
described  in  Reference  7  and  can  be  distributed  freely. 

Validation.  The  frequency  response  analysis  software  was  based  on  code 
written  for  the  Air  Force  Flight  Test  Center  (AFFTC)  CYBER  computer  system.  The 
CYBER  code  was  previously  validated  by  numerous  AFFTC  projects.  The  MATLAB™ 
version  of  the  software  was  validated  by  analyzing  data  from  several  ground  simulations 
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on  both  systems  and  comparing  the  results.  For  the  frequency  range  of  interest  (0. 1  to  10 
radians  per  second)  the  results  produced  by  the  two  versions  agreed  to  within  0. 1  percent. 

The  data  reduction  software  was  further  validated  by  running  an  example 
problem.  The  following  dynamics  were  simulated  on  MATLAB™  and  Simulink: 


_  (s  +  2)  •  e~°,25j 

Yc  (s  +  5)  (s  +  .2)  (4‘6 

A  random,  normal  input  vector  was  generated  using  the  command,  rand('normal),  and 
the  output  of  the  simulation  analyzed.  The  Bode  and  coherence  plots  are  shown  in 
Figure  C-5.  Notice  that  the  coherence  is  nearly  perfect,  especially  at  higher  frequencies. 
Additionally,  the  95%  confidence  bounds,  represented  by  the  dashed  lines  in  all  three 
plots,  are  tight.  This  frequency  response  and  the  bode  plot  of  the  actual  system  are 
shown  in  Figure  C-6.  Notice  that  except  for  the  lowest  frequency,  phase  point,  the  two 
are  nearly  identical. 


Windowing  and  Overlap.  A  common  frequency  response  analysis 
technique  is  ensemble  averaging.  Ensemble  averaging  improves  the  reliability  of  the 
frequency  response  calculations  for  time  histories  of  data  corrupted  by  measurement 
noise,  such  as  flight  test  data  (Reference  9). 

Overlap  ensemble  averaging  was  used  both  in  this  thesis  and  in 
AFFTC-TLR-93-41.  The  first  five  seconds,  or  five  hundred  data  samples,  of  each  task 
were  not  considered  so  that  the  pilot  could  become  actively  involved  in  the  task.  The 
next  4,096  samples  (40.96  seconds)  were  divided  into  seven  ensembles,  each  with  fifty 
percent  overlap.  The  data  from  several  tasks  were  analyzed  using  different  percent 
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overlap  and  numbers  of  ensembles.  In  all  cases  the  computed  frequency  response  points 
were  the  same.  The  confidence  intervals,  however,  grew  tighter  as  the  overlap  and 
number  of  ensembles  increased. 

Test  Results 

There  are  four  areas  of  analysis  that  have  implications  for  the  optimal  pilot  model 
analysis  conducted  in  this  thesis.  They  are  Pilot  Ratings  and  Comments,  Pilot  Delay, 
Statistical  Analysis,  and  Frequency  Response  Analysis.  The  following  sections  present 
the  important  results  from  each  of  these  areas. 

Pilot  Ratings  and  Comments.  The  Cooper-Harper  ratings  and  pilot  comments 
are  presented  in  Tables  C-l  and  C-2  for  the  ground  and  airborne  evaluations,  respectively. 
These  tables  also  include  the  pilot  induced  oscillation  (PIO)  rating  and  the  percentage  of 
time  the  pilot  met  desired  and  adequate  criteria  during  each  task.  Only  the  single  axis 
sum-of-sines  data  were  included  in  this  thesis. 

The  Cooper-Harper  ratings  are  presented  graphically  in  Figure  C-7.  The 
dynamics  evaluated  in  AFFTC-TLR-93-41  produced  a  wide  range  of  pilot  response,  as 
evidenced  by  the  spread  in  Cooper-Harper  ratings.  These  ratings  ranged  from  1  to  6  on 
the  ground  and  1  to  7  for  the  airborne  evaluations.  Pilot  variability  was  fairly  low.  In 
most  cases  the  Cooper-Harper  Ratings  were  within  one  rating  of  each  other.  Ratings  for 
the  airborne  evaluations  of  Case  D  were  the  exception,  ranging  from  4  to  7.  No 
significant  differences  between  ground  and  airborne  ratings  were  noted. 
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Pilot  Delay.  The  delay  between  task  command  and  pilot  response  was  estimated 
in  AFFTC-TLR-93-41  using  the  airborne,  single  axis  discrete  tracking  tasks.  This 
tracking  task  consisted  of  a  series  of  steps  and  ramps  separated  by  as  much  as  five 
seconds.  Pilot  delay  was  estimated  by  measuring  the  time  between  the  command  change 
and  the  stick  force  or  deflection  (7: 127-128). 

The  pilot  delay  between  task  command  and  stick  force  was  estimated  at  0.27 
seconds.  The  delay  times  did  not  vary  significantly  by  pilot,  dynamics  case,  or  axis. 

This  estimate  is  consistent  with  the  0.2  to  0.3  second  delay  observed  by  McRuer 
(14:217).  In  all  cases,  stick  deflection  lagged  stick  force  by  0.1  (±  0.01)  seconds, 
reflecting  the  lag  inherent  in  the  stick  dynamics.  Thus,  the  total  delay  between  task 
command  and  stick  deflection  was  0.37  seconds,  well  above  the  value  normally  used  for 
pilot  model  analysis  (7: 13-14).  [t  is  also  important  to  note  that  the  dynamics  evaluated  in 
AFFTC-TLR-93-41  were  implemented  as  a  position  command  system. 

Statistical  Analysis.  The  .andard  optimal  control  problem  minimizes  a 
performance  index  that  weights  state  deviation  and  control  usage.  The  optimal  pilot 
model  minimizes  a  performance  index  that  weights  tracking  enor  and  control  rate  usage. 
Control  rate  usage  is  included  in  the  optimal  pilot  model  performance  index  to 
implement  the  muscular  dynamics  as  described  in  Chapter  3  and  Appendix  A  of  this 
thesis  (pages  A-l,  3-5,  and  3-13).  This  section  examines  the  statistical  validity  of  both 
weightings  in  the  prediction  of  Cooper-Harper  ratings. 

Root  mean  square  (RMS)  tracking  error,  control  deflection,  and  control  rate 
values  were  computed  from  the  flight  test  data  base  produced  by  AFFTC-TLR-93-41  and 
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are  presented  in  Tables  C-3  and  C-4.  Only  the  single  axis,  sum-of-sines  data  were 
included  in  this  thesis. 

The  root  mean  square  (RMS)  values  in  Tables  C-3  and  C-4  were  computed  using 
the  following  formula: 


RMS(x)  = 


(4-7) 


In  this  equation,  n  is  the  number  of  samples  in  the  vector,  x.  The  RMS  tracking  error 
(RMS  Error)  was  determined  by  computing  the  difference  between  the  task  command 
and  the  actual  aircraft  position  (e0  =  0C  -  0)  and  then  applying  Equation  4-7  to  this  error 
vector.  The  normalized  RMS  stick  deflection  (NRMS  Stick)  was  found  by  subtracting  the 
mean  stick  deflection  from  each  element  in  the  vector  before  applying  Equation  4-7. 
Finally,  the  stick  deflection  rate  was  estimated  by  computing  the  difference  between  each 
sequential  element  in  the  stick  deflection  vector  and  dividing  the  differences  by  the 
sampling  rate,  0.01  seconds.  The  RMS  stick  deflection  rate  ( RMS  Stick  Rate)  was  found 
by  applying  Equation  4-7  to  this  stick  rate  estimate. 

The  results  of  the  regression  analysis  of  these  parameters  is  presented  in  Table  4-2 
on  the  following  page.  In  all  of  the  runs  the  first  independent  variable  was  RMS  tracking 
error  and  the  dependent  variable  was  the  Cooper-Harper  rating.  The  second  independent 
variable  was  either  normalized  RMS  stick  deflection  or  RMS  stick  deflection  rate.  The 
last  column,  variance,  was  used  to  evaluate  the  quality  of  the  least  squares  fits. 

Weighting  either  normalized  stick  deflection  or  stick  rate  produced  similar 
correlations  in  all  but  two  cases.  For  the  airborne  pitch  axis  cases,  normalized  stick 
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Table  4-2 

Regression  Analysis  of  Optimal  Pilot  Model  Weightings 


immm 


Pitch 

RMS  Error 

NRMS  Stick 

7.1936 

32.4729 

Flight 

RMS  Error 

RMS  Stick 
Rate 

4.4283 

3.9928 

Roll 

RMS  Error 

NRMS  Stick 

1.7082 

2.2877 

RMS  Error 

RMS  Stick 
Rate 

1.7371 

0.6798 

Pitch 

RMS  Error 

NRMS  Stick 

4.4451 

13.0161 

Ground 

RMS  Error 

RMS  Stick 
Rate 

1.5746 

3.4441 

Roll 

RMS  Error 

NRMS  Stick 

0.8013 

7.1414 

RMS  Error 

RMS  Stick 
Rate 

0.7280 

1.8389 

Pitch 

RMS  Error 

NRMS  Stick 

6.0179 

21.1914 

Both 
Air  and 

RMS  Error 

RMS  Stick 
Rate 

4.7125 

— . - 

3.1204 

Ground 

Roll 

RMS  Error 

NRMS  Stick 

0.9224 

7.4571 

RMS  Error 

RMS  Stick 
Rate 

1.1075 

1.4758 

Constant, 

C 

Variance, 

R2 

-12.7970 

0.8206 

-2.6624 

0.6193 

-3.9951 

0.7201 

-4.4372 

0.7233 

-3.8750 

0.5491 

0.0993 

0.9583 

-3.4535 

0.8246 

-2.7809 

0.7658 

-8.0429 

0.7311 

-2.0751 

0.7320 

-3.8638 

0.6877 

-3.8931 

0.6946 

1  The  dependent  variable  was  the  Cooper-Harper  Rating  such  that:  C-H  Rating  =  Cx,X,+  C. 


deflection  weighting  was  significantly  better  than  stick  rate  weighting.  The  reverse  was 
true  for  the  pitch  axis  ground  evaluations.  It  appears  that  both  weighting  schemes  were 
statistically  valid.  Normalized  stick  deflection  weighting  may  be  more  appropriate, 
however,  for  two  reasons.  First,  the  airborne  sample  size  was  nearly  twice  that  of  the 
ground  evaluations  (12  versus  7).  Second,  the  pilots  may  have  exhibited  singular 
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behavior  during  the  ground  simulations  due  to  the  lack  of  motion  cues.  For  example, 
they  may  have  been  more  apt  to  move  the  stick  rapidly  when  evaluating  poor  dynamics 
on  the  ground  where  there  were  no  uncomfortable  acceleration  forces  to  discourage  such 
behavior. 

The  normalized  stick  deflection  regression  coefficients  from  the  airborne 
evaluations  were  used  to  draw  the  regression  analysis  lines  in  Figures  C-8  and  C-9.  Each 
airborne  test  point  was  plotted  as  a  function  of  RMS  error  and  normalized  RMS  stick 
deflection.  The  number  above  each  data  point  is  the  Cooper-Harper  rating  assigned 
during  the  task.  As  shown  in  these  figures,  there  was  a  strong  relationship  between  these 
two  variables  and  the  Cooper-Harper  rating,  but  this  relationship  was  by  no  means 
perfect. 

Frequency  Response  Analysis 

From  a  linear  systems  perspective,  the  primary  input  to  the  pilot  was  assumed  to 
be  task  error.  The  pilot's  primary  outputs  were  assumed  to  be  stick  force  and  deflection. 
A  frequency  response  analysis  of  these  input  and  output  variables  was  performed  for  all 
of  the  single  axis,  sum-of-sines  tracking  tasks.  This  section  presents  the  applicable 
results  of  this  analysis. 

Stick  Force  Versus  Stick  Deflection.  As  described  previously  in  this 
chapter,  the  evaluated  dynamics  were  implemented  as  a  position  command  system.  Stick 
force  and  displacement  were  related  by  a  high  frequency  stick  dynamics  term  and  a  force 
gradient  as  shown  in  Figure  C-l .  The  only  difference  between  the  frequency  responses 
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of  stick  force  and  stick  deflection  was  the  frequency  response  of  these  stick  dynamics.  In 
other  words,  when  the  frequency  responses  of  stick  force  to  task  error  were  summed  with 
the  Bode  plot  of  the  stick  dynamics,  the  result  was  identical  to  the  stick  deflection 
frequency  responses  in  Appendix  D.  As  a  result,  the  same  insights  could  have  been 
gained  by  studying  either  variable.  The  analysis  in  AFFTC-TLR-93-41  and  this  thesis 
focused  on  stick  deflection  (7: 14-15). 

Ground  Simulation.  As  discussed  in  AFFTC-TLR-93-41,  there  were  no 
significant  differences  between  the  pilot's  airborne  and  ground  frequency  response.  The 
phase  lead  at  higher  frequencies,  however,  was  slightly  greater  on  the  ground  than  in  the 
air.  There  were  no  consistent  trends  for  the  magnitude  (7:1 8). 

Stick  Deflection  to  Task  Error,  A  frequency  response  analysis  of  stick 
deflection  to  tracking  task  error  for  all  single  axis  sum-of-sines  tracking  tasks  was 
completed  for  AFFTC-TLR-93-41 .  A  representative  cross-section  of  this  data  is 
presented  in  Appendix  D. 

The  power  spectral  densities  (PSD)  of  the  pilot  input  and  the  pilot  output  for  two 
of  the  worst  dynamics  cases,  Cases  4  and  D,  are  presented  in  Figures  D-l  and  D-2.  Even 
for  these  two  highly  oscillatory  cases  there  were  no  notches  in  the  PSD  magnitude. 

Likewise,  the  Bode  plots  of  the  pilot  response  did  not  reveal  any  higher  order 
compensation,  or  notching  behavior.  One  representative  frequency  response  plot  for  each 
dynamics  case  is  presented  in  Figures  D-3  through  D-10.  As  shown  in  these  figures,  the 
magnitudes  were  consistent  with  the  classical  gain,  lead,  and  lag  pilot  model  form 
described  previously  in  this  thesis  and  in  Reference  15.  The  phase,  however,  was  not 
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consistent  with  this  form  of  pilot  compensation  due  to  the  large  amount  of  phase  lead 
present  at  higher  frequencies.  As  shown  in  Figure  D-10,  this  high  frequency  phase  lead 
is  especially  noticeable  for  cases  with  delay.  Despite  this  phase  lead,  no  higher  order 
compensation  was  noted  in  any  of  the  evaluations. 

The  pilot  response  appears  to  be  strongly  related  to  the  aircraft  dynamics.  To 
see  this  relationship,  the  frequency  responses  of  the  combined  pilot-aircraft  systems  are 
presented  in  Figures  D-l  1  through  D-1S.  These  plots  were  formed  by  adding  the 
frequency  response  of  the  pilot  (stick  deflection  to  task  error)  with  that  of  the  aircraft 
(pitch  or  roll  angle  to  stick  deflection).  As  shown  in  these  figures,  the  combined  system 
responses  resembled  integrators  near  the  cross-over  frequency  in  all  cases.  This  is 
consistent  with  McRuer's  crossover  pilot  model  theory  (Reference  17).  Note  also  that  the 
cross-over  frequency  of  the  combined  system  dropped  from  1.9  radians  per  second  for 
the  best  pitch  case  (Figure  D-l  1)  to  1.6  radians  per  second  for  the  worst  pitch  case 
(Figure  D-14)  and  from  1.6  radians  per  second  for  the  best  roll  case  (Figure  D-l 5)  to 
1.0  radians  per  second  for  the  worst  roll  case  (Figure  D-l 8). 

Gain  Effects.  The  effect  of  command  path  gain  was  briefly  evaluated 
during  the  ground  simulation  phase.  The  frequency  response  analysis  of  stick  deflection 
to  task  error  for  the  ground  evaluation  of  Case  1  is  shown  in  Figure  D-l 9.  The  command 
path  gain  for  Case  1  was  increased  from  8  to  1 6  and  re-evaluated  by  the  same  pilot.  The 
resulting  frequency  response  is  shown  in  Figure  D-20.  As  shown  in  this  figure,  the  pilot 
decreased  his  gain  by  about  6  dB  to  compensate.  This  behavior  is  even  more  apparent 
when  examining  the  combined  pilot-aircraft  systems  in  Figures  D-21  and  D-22.  Even 
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though  the  aircraft  command  path  gain  for  Figure  D-22  was  twice  that  in  Figure  D-21, 
the  frequency  responses  of  the  combined  systems  were  nearly  the  same. 

Conventional  Pilot  Model  Predictions 

The  dynamics  cases  described  in  Appendix  C  were  evaluated  using  the  STI 
optimal  pilot  model  and  all  of  the  applicable  pitch  and  roll  axis  pilot  models  in 
MIL-STD-1797  (Reference  5). 

The  results  of  the  pitch  axis  pilot  model  analysis  are  presented  in  Table  C-6.  As 
shown  in  this  table,  the  MIL-STD-1797A  pitch  models  were  unsatisfactory  for  predicting 
the  handling  qualities  ratings  of  the  pitch  dynamics  used  for  this  evaluation.  Three  of  the 
five  models  predicted  Level  II  handling  qualities  for  Case  1  which  actually  had  Level  I 
handling  qualities.  All  of  the  models  predicted  Level  III  handling  qualities  for  Cases  3 
and  4,  and  both  of  these  cases  had  Level  II  handling  qualities. 

The  optimal  pilot  model  was  moderately  successful  in  predicting  the 
Cooper-Harper  ratings,  but  the  predicted  pilot  describing  functions  were  not  consistent 
with  the  observed  behavior.  Bode  plots  of  the  predicted  pilot  describing  functions  and 
the  resulting  pilot-aircraft  systems  are  shown  in  Figures  C-10  through  C-13  for  Cases  1 
and  4.  Due  to  the  higher  order  pilot  compensation,  or  notching,  present  in  these  figures, 
they  do  not  resemble  the  corresponding  flight  test  responses  in  Figures  D-3,  D-6,  D-ll, 
and  D-14. 

The  roll  axis  pilot  model  predictions  are  presented  in  Table  C-7.  These  models 
adequately  predicted  the  Cooper-Harper  ratings  and  handling  qualities  levels  of  the  roll 
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axis  dynamics.  Due  to  the  simplicity  of  the  MIL-STD-1797A  roll  axis  models,  however, 
these  models  had  to  be  used  together  in  a  checklist  fashion  to  gain  adequate  insight  into 
the  aircraft's  predicted  handling  qualities.  The  optimal  pilot  model  predictions  were 
slightly  pessimistic,  while  the  bandwidth  criterion  predictions  were  slightly  optimistic. 
Both  were  accurate  enough  to  be  useful.  As  with  the  pitch  axis,  the  pilot  describing 
functions  predicted  by  the  STI  optimal  pilot  model  were  not  consistent  with  observed 
behavior. 

Summary 

In  an  effort  to  provide  insight  into  human  pilot  behavior,  a  limited  evaluation  of 
human  pilot  response  was  conducted.  This  evaluation  consisted  of  ground  and  airborne 
simulations  in  the  variable  stability  Lear  II  aircraft.  Four  different  pitch  and  four 
different  roll  axis  dynamics  were  evaluated  using  three  different  tracking  tasks.  For  each 
set  of  dynamics,  primary  pilot  response  parameters  were  recorded  and  examined  using 
Fourier  transform  analysis. 

The  dynamics  evaluated  produced  a  wide  range  of  pilot  response  as  evidenced  by 
the  spread  in  Cooper-Harper  ratings.  Pilot  variability  was  fairly  low.  In  all  but  one  case 
the  Cooper-Harper  ratings  were  within  one  rating  of  each  other.  No  significant 
difference  between  ground  and  airborne  ratings  were  noted  (page  4-12)1. 

Using  the  discrete  tracking  task  time  histories,  the  pilot  delay  between  task 
command  and  stick  force  was  estimated  at  0.27  seconds.  In  all  cases,  stick  displacement 
lagged  stick  force  by  0. 1  seconds  due  to  the  lag  inherent  in  the  stick  dynamics.  Thus,  the 
1  Refers  to  the  page  in  this  thesis  containing  this  conclusion. 
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total  delay  between  command  and  stick  displacement  was  0.37  seconds,  well  above  the 
value  normally  used  for  pilot  model  analysis  (page  4-13). 

The  statistical  validity  of  using  root  mean  square  (RMS)  tracking  error  and  stick 
deflection  or  deflection  rate  to  predict  Cooper-Harper  ratings  was  examined.  A 
regression  analysis  revealed  that  task  error  and  normalized  stick  deflection  weighting 
produced  the  best  correlation  with  airborne  pilot  ratings  (page  4-15). 

A  frequency  response  analysis  of  stick  force  and  displacement  to  task  error  was 
conducted  for  all  single  axis,  sum-of-sines  tracking  tasks.  The  following  conclusions 
were  made. 

1 .  The  only  difference  between  the  frequency  responses  of  stick  force  and  stick 
deflection  was  the  frequency  response  of  the  stick  dynamics.  Thus,  the  same 
insights  could  have  been  gained  by  examining  either  variable.  Because  the 
evaluated  dynamics  were  implemented  as  a  position  command  system,  stick 
deflection  was  chosen  (page  4-16). 

2.  There  were  no  significant  differences  between  the  pilot's  airborne  and  ground 
frequency  response  (page  4-17). 

3.  No  notches  were  observed  in  any  of  the  stick  displacement  power  spectral 
densities  (page  4-17). 

4.  No  higher  order  pilot  compensation,  or  notching  behavior,  was  noted  in  any 
of  the  single  axis  frequency  responses.  The  magnitude  plots  in  both  axes 
were  consistent  with  the  classical  gain,  lead,  and  lag  pilot  model  form.  The 
phases  were  not  consistent  with  this  form  due  to  large  amounts  of  pure  phase 
lead  present  at  higher  frequencies  (page  4-17). 

5.  In  all  cases,  the  pilot  applied  compensation  so  that  the  combined  pilot-aircraft 
system  resembled  an  integrator  near  the  cross-over  frequency  (page  4-18). 

6.  When  the  aircraft  command  path  gain  was  doubled,  the  pilot  reduced  his  gain 
so  that  the  response  of  the  combined  pilot-aircraft  system  remained 
unchanged  (page  4-18). 
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All  of  the  dynamics  were  evaluated  by  the  pilot  models  in  MIL-STD-1797A  and 
by  the  STI  optimal  pilot  model.  The  MIL-STD-1797A  roll  axis  pilot  models  adequately 
predicted  the  Cooper-Harper  ratings  and  handling  qualities  levels  of  the  dynamic  used  for 
this  project.  The  MIL-STD-1797A  pitch  axis  models,  however,  were  unsatisfactory  for 
predicting  the  handling  qualities  levels  or  ratings.  The  STI  optimal  pilot  model 
successfully  predicted  the  Cooper-Harper  ratings  in  both  axes.  The  pilot  describing 
functions  predicted  by  this  model  were  not  consistent  with  observed  behavior  because 
they  contained  higher  order  compensation  (page  4-19). 
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5.  A  Numerical  Solution  to  the  Linear  Quadratic  Gaussian  Problem 
General 

The  standard  solution  to  the  linear  quadratic  Gaussian  (LQG)  problem  is  a 
compensator  of  the  same  order  as  the  controlled  element.  Because  the  optimal  pilot 
model  is  based  on  LQG  theory,  the  high  order  pilot  describing  functions  it  predicts  .ire 
generally  not  consistent  with  classical  pilot  modeling  theory  or  observed  behavior.  The 
method  detailed  in  this  chapter  was  developed  to  provide  a  way  around  this  problem. 

The  numerical  LQG  solution  described  in  this  chapter  allows  the  user  to  select  a 
compensator  form.  It  then  uses  a  Nelder-Meade  simplex  algorithm  to  find  the 
coefficients  of  the  compensator  that  minimize  the  standard  LQG  performance  index. 

This  method  is  not  only  useful  in  the  optimal  pilot  model  problem.  It  may  apply  to  all 
situations  when  reduced  order  compensation  is  desired. 

This  chapter  is  divided  into  four  section.  The  first  section  of  this  chapter  gives 
general  background  information.  The  second  section  describes  the  numerical  LQG 
solution  method.  This  method  is  then  demonstrated  by  solving  two  example  problems. 
Finally,  the  results  of  this  chapter  are  summarized. 

Background 

Reduced  Order  Controllers.  There  are  three  principal  approaches  to  reduced 
order  LQG  controller  design  as  shown  in  Figure  5-1  (1:253).  One  method  is  to  reduce 
the  order  of  the  controlled  element.  For  the  optimal  pilot  model  problem  this  alone  is  not 
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Figure  5-1.  Reduced  Order  Controller  Design 


a  practical  solution.  The  second  approach  requires  the  application  of  a  controller 
reduction  algorithm  to  the  standard  LQG  solution.  As  shown  in  Chapter  3  and  Appendix 
B,  however,  accurate  first  order  approximations  of  optimal  pilot  model  results  are  not 
always  attainable.  The  more  direct  approach  is  to  use  a  constrained  optimal  control 
methodology  similar  to  the  one  described  in  this  chapter.  For  a  detailed  discussion  of 
constrained  optimization  theory  consult  Bernstein  and  Hyland  (Reference  1). 

The  Standard  Linear  Quadratic  Gaussian  Problem .  A  general  form  of  the 


standard  LQG  problem  is  shown  in  Figure  5-2. 


Figure  5-2.  General  Linear  Quadratic  Gaussian  Problem 
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where 


Bc,  Cc 

Controlled  Element  State  Space  Matrices 

xc 

Controlled  Element  State  Vector 

AK <  B* ,  Ljj,  Dk 

Compensator  State  Space  Matrices 

xK 

Compensator  State  Vector 

u 

Control  Vector 

y 

Controlled  Element  Output  Vector 

ri(t) 

Measurement  Noise 

m 

Disturbance  Noise 

r(s) 

Disturbance  Noise  Coloring  Filter  (Laplace) 

The  methods  developed  in  this  thesis  assumed  the  controlled  element  was  time-invariant 
and  strictly  proper  (no  Dc  term)1.  These  assumptions  were  not  considered  limiting. 

The  standard  LQG  solution  is  the  compensator  that  minimizes  the  following 
performance  index: 

ao 

J=  jlxTQx  +  uTRu\dt  (5-1) 

o 

Where  R  is  the  control  weighting  and  Q  is  the  state  deviation  weighting.  If  the  system  is 
driven  only  by  random  noises  as  shown  in  Figure  5-2,  this  performance  index  can  be 
rewritten  as  an  expected  value  rather  than  an  integral  as  shown. 


J  =  E{xTQx  +  uTRu )  (5-2) 

From  Figure  5-2,  the  control  vector,  u,  is: 


u  -  -[Ck  xk  +  Dk  CcXc  +  Dkx\\ 


(5-3) 


Notice  that  the  last  part  of  this  expression,  DKx\,  contains  a  linear  combination  of  the 
sensor  noises.  If  this  expression  is  substituted  into  Equation  5-2,  the  performance  index 


1  This  assumption  was  made  to  simplify  the  numerical  solution  derivation. 
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will  contain  the  following  term: 


E{x\TDTKRDKr\}  (5-4) 

Unless  this  value  is  zero,  the  performance  index  will  be  infinite.  For  illustration 
purposes,  assume  there  is  only  one  feedback  channel.  The  matrices,  DK  and  R,  become 
scalars,  and  Equation  5-4  can  be  written  as: 

D\  R  E{ n2}  (5-5) 

The  sensor  noises  in  the  LQG  problem  are  defined  as: 

E{ri(r)-Tir(/-t»  =R0-5(x)  Vf,t  (5-6) 

where  R„  is  the  intensity  of  the  measurement  noise.  Combining  terms,  Equation  5-5  can 
be  written  as: 

D2kRR0-  5(t)  (5-7) 

Because  of  the  delta  function,  the  magnitude  of  Equation  5-7  is  infinite.  Thus,  the 
magnitude  of  the  performance  index,  which  contains  this  term,  is  also  infinite. 

The  most  obvious  solution  to  this  problem  is  to  require  the  compensator  feed-through 
term,  DK,  to  be  zero.  In  other  words,  the  compensator  must  be  strictly  proper. 

Assuming  a  proper  compensator,  the  LQG  problem  can  be  represented  by  the 
block  diagram  in  Figure  5-3.  Also,  the  frequency  content  of  the  disturbance  noise,  T(s), 
is  now  written  in  state  space  form. 
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Figure  5-3.  Linear  Quadratic  Gaussian  Problem 
—  Proper  Compensator 


Assuming  the  controlled  element  disturbance  noise  filter  is  also  strictly  proper',  the 
equations  of  state  for  this  system  are: 


Xc  =AcXc+Bctt  +  CrXr 

y  =  CcXc 

(5-8) 

Xr  =AtXy  +2?r^ 

(5-9) 

xK  =  AKXK+BK(r\  +7) 

u  =  -Ck  xk 

(5-10) 

Substituting  and  combining  terms  yields  the  following  state-space  representation  driven 
only  by  white  noises: 


xc 

Ac 

Cr  -BcCk 

xc 

*  0 

0 

.  n  . 

xr 

.  . 

— 

0 

.  BkCc 

Ar  0 

0  Ak 

xr 

.  **  . 

+ 

B  r 
0 

£  ° 

1 _ 

1  This  assumption  was  made  to  simplify  the  numerical  solution  derivation. 


(5-11) 
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This  system  is  in  the  form  of: 


x  =  Ax  +  Ew 


where  the  noise  intensities  are: 

£{5(o-5r(»-t)}=e.-6(x) 
£(1(0  ■  nr(<-  ')}=£.  ■  S(T) 
and  the  constant  intensity  noise  matrix,  Q„,  is  defined  as: 


(5-12) 


(5-13) 

(5-14) 


(5-15) 


One  important  point  should  be  made.  When  there  is  no  driving  noise  coloring 
filter,  the  filter  states  should  be  removed.  If  Ar ,  Br ,  and  Cr  are  simply  set  to  zero,  the 
matrices  in  Equation  5-11  will  not  form  a  minimum  realization.  The  state  space 
representation  of  the  system  without  coloring  filters  is  given  in  Equation  5-16. 


xc 

XK 


Ac  -BcCk 
BkCc  Ak 


xc 

xk 


+ 


r  o 

0  Bk 


(5-16) 


Where  T  is  now  a  gain  matrix  that  distributes  the  driving  noise  into  the  controlled 
element  states. 


The  Numerical  Solution 

As  mentioned  earlier,  for  systems  driven  by  random  noises,  the  standard  LQG 
performance  index  can  be  written  as  the  sum  of  two  expected  values. 

J=E{xTQx)  +E{uTRu }  (5-17) 
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The  expected  value  of  the  state  vector  can  be  found  from  the  covariance  matrix,  X,  such 
that: 

X=E{xxT)  (5-18) 

where  the  covariance  matrix  is  the  solution  to  the  following  Lyapunov  equation 
(11:106-111). 

AX+XAt+EQ„Et  =  0  (5-19) 

The  matrices,  E  and  A,  correspond  to  those  in  Equation  5-12,  and  Qn  is  the  noise  matrix 
defined  in  Equation  5-15. 

The  covariance  matrix  can  be  used  to  find  the  value  of  the  performance  index. 
Equation  5-17,  by  defining  the  following  matrices. 


XG  =  [In  \0nxm]X=  ' C 

(5-20) 

XfC  —  [Oirijtn  1  Im  ] X 

(5-21) 

U  =  -CrXk 

(5-22) 

where 

xG=  Forward  Loop  State  Vector  (Controlled  Element  and  Filter) 
n  =  Number  of  Controlled  Element  and  Filter  States 
m  =  Number  of  Compensator  States 
I  =  Identity  Matrix 

Combining  Equations  5-20  to  5-22  and  Equation  5-17  yields: 

J=E{xtg  Qxg }  +  E{xTK  CtkR  Ckxk }  (5-23) 
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Assuming  Q  is  a  diagonal  matrix,  the  first  term  of  Equation  5-23  can  be  expanded  as 
follows. 


E{xTqQXg)  =  E{XgiQ\\Xg\  +XG2Q22XG2  +  ..JCGiQtiXGi)  (5-24) 


Grouping  terms,  Equation  5-23  can  be  written  as: 


E{xIQxg]  =  X  Qu-E{x2Gi} 

i*  1 


=  tQirXu  (5-25) 

<=i 


where 

n  =  Number  of  Controlled  Element  and  Filter  States 
Q  =  State  Deviation  Weighting  Matrix  from  Equation  5-17 
X  =  Covariance  Matrix  from  Equation  5-18 

The  second  term  of  Equation  5-23  is  a  little  more  convoluted.  If,  for  example, 
there  are  two  compensator  states  and  two  feedback  channels,  this  term  expands  to: 


Xki 

CkU  Ck21 

^11 

0 

Cku 

Ck\2 

xki 

. Xfc2 

.  CkU  Ck22 

0 

R-22 

.  C/C21 

CfC22 

.  Xki 

=  Eh\_XkiC2ku  +xjr2  C2Kx2  +2xkiXkiCkuCku~\  + 


Ez^x2KXC2K2X  +Xk2C2K22  +  2xkiXk2Ck2iCk22]  (5-26) 


Grouping  terms,  the  general  form  of  this  expression  is: 


E{xTkCTkRCkXk]  =  X  X  X  Rhh-  CK(ht)'CK(hf)  ■  E{XK(l}  XK(f)} 

(=■=  1  j=  1 


-  X  X  X  Rfik  ‘  Cfc(hf)  '  Cx(hf)  ‘  -^(ii +<,»+/)  (5-27) 

fc=l  i=lj=\ 
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where 


o  =  Length  of  Control  Vector,  u  (number  of  rows  in  CK) 
m  =  Number  of  Compensator  States 
n  =  Number  of  Controlled  Element  and  Filter  States 
R  =  Control  Weighting  from  Equation  5-17  (Diagonal) 

CK=  Compensator  State  Space  Matrix 
X  =  Covariance  Matrix  from  Equation  5-18 

Equations  5-25  and  5-27  are  the  central  equations  in  the  numerical  solution  of  the 
LQG  problem.  Using  these  equations  the  performance  index  can  be  computed  for  any 
controlled  element  and  compensator.  A  numerical  minimization  routine  can  then  search 
all  of  the  compensator  parameters  for  the  combination  that  returns  the  minimum 
performance  index.  Additionally,  the  performance  index  can  be  modified  to  include  side 
constraints,  such  as  gain  or  phase  margin. 

Equations  5-25  and  5-27  are  implemented  in  the  MATLAB™  routines,  or  ".M" 
files  listed  in  Appendix  E.  PIFIND.M  computes  the  performance  index  value  for 
problems  with  driving  noise  coloring  filters.  PIFINDNF.M  computes  the  performance 
index  value  for  problems  without  coloring  filters. 

Example  Problems 

This  numerical  solution  is  illustrated  in  the  following  examples.  The  first 
example  consists  of  a  low  order  plant  augmented  with  a  disturbance  noise  coloring  filter. 
The  second  example  uses  a  high  order  plant  with  no  disturbance  noise  coloring  filter. 
These  examples  are  not  necessarily  meant  to  mimic  aircraft  dynamics.  This  numerical 
method  can  apply  to  any  type  of  LQG  problem. 
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Low  Order  Plant  with  Coloring  Filter.  Assume  the  controlled  element  is  an 


integrator  and  the  driving  noise  is  filtered  as  follows: 

State  Space  Representation 

Yc  =  j  =>  [xc]  =  [0]xr  +  [l]w 

y  =  [\\xc  (5-28) 

Fr  =  y  =  -jU  =>  {*r]  =  [-2]xr  +  UH 

W'c  =  [  TO”  ]  Xr  (5-29) 


where  %  is  the  random  disturbance  noise  and  Wc  is  this  noise  after  the  coloring  filter. 
Assume  the  state  deviation  weighting  is  one  and  the  control  usage  weighting  is  two.  The 
state  deviation  weighting  matrix,  Q,  must  be  diagonal  due  to  the  assumptions  made  in  the 
development  of  the  numerical  solution.  Because  the  controlled  element  states  must  be 
augmented  with  the  filter  states,  the  state  deviation  weighting  matrix,  Q,  must  have  the 
same  number  of  rows  as  the  total  number  of  filter  and  controlled  element  states.  There  is 
normally  no  penalty  on  filter  state  deviations,  making  their  weightings  zero.  The  control 
usage  weighting  matrix,  R,  must  be  diagonal  with  the  same  number  of  rows  as  in  the 
control  vector,  u.  In  this  example,  there  are  two  total  controlled  element  and  filter  states 
and  only  one  control  channel.  The  weighting  matrices  are  shown  in  Equation  5-30. 


Q  = 


1 

0 


0 

0 


*  =  [>] 


(5-30) 
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The  noises,  t,  and  r\,  will  have  unit  intensity  for  this  example  such  that  the  noise  intensity 
matrix,  QH,  is  as  shown  in  Equation  5-31. 


Qn  = 


1  0 
0  1 


(5-31) 


In  this  example,  assume  the  optimal  second  order  compensator  of  the  form  in 
Equation  5-32  is  desired. 


K()  _  c(l)  s  +  c(2) 
s2  +  c(3)  •  s  +  c(4) 


(5-32) 


The  Nelder-Meade  simplex  algorithm,  finins,  on  MATLAB™  will  now  be  used  to  search 
for  the  combination  of  compensator  variables,  c(i),  that  produces  the  minimum 
performance  index.  An  additional  function  (. m  file )  containing  the  compensator  form 
must  be  written.  The  routine  for  solving  this  problem  is: 

fimction[J]=comp(c,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise) 

[Ak,Bk,Ck,Dk]=tf2ss([c(l),c(2)],abs([l,c(3),c(4)])); 

J=PIFIND(Ak,Bk,Ck,Ac,Bc,Cc,Ag,Bg»Cg,Q,R,Noise); 

The  first  line  lists  the  input  and  output  arguments.  The  second  line  converts  the  desired 
compensator  form  into  state  space,  and  the  final  line  finds  the  performance  index,  J, 
using  the  file,  PIFIND,  listed  in  Appendix  E.  This  function,  comp,  is  minimized  by 
typing  the  following  MATLAB™  command. 

cmin=fminsCcomp',cO,[3,Q,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise) 

In  this  command,  cO  is  the  initial  guess  and  the  rest  of  the  matrices  are  as  shown  in 
Equations  5-28  through  5-3 1 .  Noise  refers  to  the  matrix,  Qn ,  shown  in  Equation  5-3 1 . 
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Due  to  the  nature  of  the  Nelder-Meade  simplex  search  routine,  the  initial  guess  requires 
some  consideration.  If  any  evaluated  compensator  has  poles  in  the  right  half  plane  or  on 
the  imaginary  axis,  the  simplex  search  may  not  be  successful.  Thus,  restricting  the 
denominator  values  in  the  search  vector,  c,  to  positive  values  greatly  improves  the 
robustness  of  the  numerical  search.  Likewise,  the  Lyapunov  solution  must  be  unique. 

The  necessary  and  sufficient  condition  for  this  is: 

+  \j(A)  *  0  for  ViJ=l...n  (5-33) 

where  A  is  from  Equation  5-12  and  A  is  the  eigenvector  of  A.  A  sufficient  condition  for  a 
unique  Lyapunov  solution  is  that  the  real  part  of  all  eigenvectors  be  negative  (3:30). 

Thus,  the  search  may  not  converge  for  compensators,  controlled  elements,  or  filter 
dynamics  with  entire  state  space  matrices  of  zeros.  As  the  number  of  parameters  to 
search  increases,  the  importance  of  the  initial  guess  increases.  Additionally,  the 
Nelder-Meade  routine  is  best  for  finding  the  minimum  of  functions  with  five  or  fewer 
unknown  parameters  (4:1 16-122).  For  complex  compensator  forms  a  different  numerical 
search  routine  may  be  appropriate.  For  this  example,  any  initial  condition  vector,  cO, 
with  values  between  approximately  0.1  and  100  will  return  the  proper  solution. 

In  this  example,  the  initial  guess  was 

YKi  =  —  ~  or  cO  =  [  1  1  1  1  ]  (5-34) 

s2+s+l  L  J 

and  the  following  compensator  was  computed: 

Cmin  =[  0.9876  2.0977  3.8588  4.2075  ]  (5-35) 
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(5-36) 


Y  0.9876^  +  2.0977 
*  s2  +  3.85885  +  4.2075 

with  a  resulting  performance  index,  J,  of  3.9252. 

As  a  check,  the  same  problem  was  worked  using  the  conventional  LQG  solution 

methods  described  in  Appendix  A.  First,  the  system  had  to  be  put  into  the  following 

form. 


x  =  Ax  +  Bu  +  Tt, 

y  =  Cx+Dr\  (5-37) 

This  was  accomplished  using  the  following  formulas. 


Xc 
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■u  + 

'  0 

. ir  - 

.  o  At  . 
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T 
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_Br  . 

(5-38) 


r=[Cc  0]. 


Xc 

xr 


+  [0]«  +  r| 


(5-39) 


For  this  example  these  equations  were: 


xc 

xr 


0  JiTs 
0  -2 


xc 


■u  + 


0 

1 


(5-40) 


Y=[  1  0]- 


xc 

xr 


+  [0]  •  u  +  r\ 


(5-41) 


In  the  previous  method,  the  system  was  only  driven  by  random  noises  and  was  in  the 
form  of  Equation  5-12.  Now  the  system  is  driven  both  by  noise  and  the  control  input,  u. 
Because  of  this,  the  disturbance  noise  matrix  is  now: 


Qf=r-Q0  rT 


(5-42) 


5-13 


where  Q0  is  the  noise  intensity  from  Equation  5-13,  and  T  is  the  disturbance  distribution 
gain  matrix  in  Equation  5-37.  The  state  deviation  matrix,  Q,  control  weighting  matrix,  R, 
and  measurement  noise  matrix,  R0 ,  are  the  same  as  in  the  numerical  solution. 

For  the  Iqg  command  in  MATL  AB™,  the  control  weightings  and  noises  must  be 
grouped  as  shown  in  Equations  5-43  and  5-44. 


W= 


Q  0 

0  R 


1  0  0 
0  0  0 
0  0  2 


V= 


Qf  o 

0  R0 


0  0  0 
0  1  0 
0  0  1 


(5-43) 


(5-44) 


The  following  command  returned  the  same  solution  as  that  found  in  Equation  5-36  using 
the  numerical  method. 

[Ak,Bk,Ck,Dk]=lqg(A,B,C,D,W,V) 

where  A,  B,  C,  and  D  are  from  Equations  5-40  and  5-41  and  fV and  Fare  from  Equations 
5-43  and  5-44. 

High  Order  Plant  with  No  Coloring  Filter.  This  example  is  presented  to 
illustrate  the  application  of  this  method  to  reduced  order  compensator  design.  In  this 
example  the  following  controlled  element  state  space  was  used: 
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"  1  " 
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(5-45) 
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y  =  [  0  1  6  8  ]  •  xc  +  [0]  •  u 


(5-46) 


This  is  a  single-input-single-output  system  (SISO)  as  shown.  Unit  intensity  weightings 
were  applied  to  the  control  usage  and  state  deviations  such  that: 


R  =  [  1  ]  and  Q  = 


10  0  0 
0  10  0 
0  0  10 
0  0  0  1 


(5-47) 


The  control  and  state  deviation  weighting  matrices,  R  and  Q,  must  be  diagonal  for  the 
numerical  solution  to  be  valid.  The  noises,  E,  and  x\,  were  assumed  to  be  unit  intensity, 
random  white  noises  such  that: 


- 1 

to 

o 

O 

_ i 

’  1  0  ' 

l  0  Ro 

0  1 

(5-48) 


and  the  standard  LQG  noise  intensity  matrix  was: 


Qf=TQ0TT  = 


10  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


(5-49) 


Using  the  MATLAB™  command,  Iqg,  described  in  the  previous  example,  the  resulting 


standard  LQG  solution  was: 


YkLeG(s)  = 


-0.0163s3 -0.1413  s2 -0.19295- .0680 
s4  +4.7837  s3  + 12.9420  s2  +  17.8484  s  + 9.7295 


(4-50) 
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The  numerical  LQG  solutions  were  found  using  the  same  method  as  in  the 
previous  example.  The  routine  PIFINDNF,  listed  in  Appendix  E,  computed  the 
compensators  listed  in  Table  5- 1 . 


Table  5-1 

Sub-Optimal  Compensator  Example 


Order 

Compensator 

PI1 

Time2,3 

4 

-0.0163s3  -  0. 141 5  s2  -  0. 1949  s  -  .0689 
s4 +  4.7981  s3  + 13.0013  s2 +  17.9841  s  + 9.8708 

0.23694 

2:20 

3 

-0.0162s2 -0.1284s -0.0620 
s3  +  3.8922  s2  +  9.0063  s  +  8.9584 

0.23694 

1:57 

2 

-0.0237s -0.0398 
s2  + 1.3466s +  4.3805 

0.23695 

0:45 

1 

-0.0335 
s  + 2.6522 

0.23713 

0:10 

1  Performance  Index 

2  Processing  time  (minutes:  seconds)  on  IBM  compatible  486DX  (50  MHz) 

3  Computing  time  for  the  standard  LQG  solution  method  was  five  seconds  on  the  same  computer. 


Two  conclusions  can  be  drawn  from  these  examples.  First,  the  performance  index 
computation  routines  are  valid.  When  the  desired  compensator  order  was  equal  to  the 
LQG  compensator  order,  the  numerical  method  and  the  standard  LQG  solution  method 
produced  the  same  results,  within  the  tolerance  of  the  numerical  search.  Second, 
reducing  the  compensator  order  did  not  significantly  affect  the  resulting  performance 
index  value.  For  less  than  one  tenth  of  a  percent  increase  in  the  performance  index,  a 
first  order  compensator  could  have  been  used  in  this  example  instead  of  the  fourth  order, 
standard  LQG  solution. 
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For  higher  order  compensators  (more  than  six  unknown  variables),  this  simplex 
search  algorithm  was  not  satisfactory  and  an  alternative  numerical  minimization  routine 
may  be  necessary.  As  the  number  of  unknown  variables  increased  above  six,  the  solution 
became  sensitive  to  the  initial  guess  for  some  problems.  For  the  second  example 
problem,  the  third  order  solution  was  found  using  the  following  initial  condition: 

Co  =  [  -0.1  -0.1  -0.1  4  4  4  ]  (5-51) 

The  fourth  order  solution  was  found  using  a  small  perturbation  of  the  standard  LQG 
solution.  If  a  higher  order  compensator  is  desired,  the  MATLAB™  Lyapunov  routine 
should  be  rewritten  so  that  it  does  not  return  an  error  message  and  interrupt  the  program 
when  the  Lyapunov  solution  is  not  unique.  A  large,  arbitrary  performance  index  should 
be  assigned  instead.  The  Nelder-Meade  simplex  algorithm  was  satisfactory  for  the  lower 
order  compensators.  The  first  and  second  order  compensators  were  not  sensitive  to  the 
initial  guesses.  Any  initial  condition  with  values  between  1  and  WO  allowed 
convergence  to  the  same  optimal  solution.  Additionally,  the  computing  time  for  the  first 
and  second  order  compensators  was  not  significant. 

Bode  plots  of  the  four  compensators  are  shown  in  Figures  E-l  through  E-10.  As 
shown  in  Figures  E-9  and  E-10,  the  second,  third,  and  fourth  order  compensators  are 
roughly  equivalent.  From  a  low  order  equivalent  system  standpoint,  the  first  order 
compensator  was  not  a  good  approximation  of  the  others.  However,  this  compensator 
only  increased  the  LQG  performance  index  0.08  percent  over  the  fourth  order  solution  as 
shown  in  Table  5-1. 
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Summary 


A  numerical  solution  to  the  LQG  problem  was  developed.  This  method  allowed 
the  user  to  select  a  compensator  form.  It  then  used  a  Nelder-Meade  simplex  algorithm  to 
find  the  coefficients  of  the  compensator  that  minimized  the  standard  LQG  performance 
index.  This  method  is  valid  for  situations  when  reduced  order  LQG  compensation  is 
desired. 

Several  assumptions  were  made  in  the  development  of  the  numerical  LQG 
solution. 

1 .  The  controlled  element  is  time  invariant  (page  5-3)1. 

2.  The  controlled  element  is  proper  (page  5-3). 

3.  The  compensator  is  proper  (page  5-4). 

4.  The  system  is  driven  only  by  random  noises  (page  5-3). 

5.  The  frequency  content  of  the  driving  noise  is  proper  (page  5-5). 

6.  The  state  deviation  weighting  matrix,  Q,  is  diagonal  (page  5-8). 

7.  The  control  usage  weighting  matrix,  R,  is  diagonal  (page  5-9). 

Given  these  assumptions.  Equations  5-25  and  5-27  were  derived.  These  equations 
compute  the  LQG  performance  index  value  for  any  controlled  element  and  compensator. 
Two  routines  that  implement  these  equations  were  provided  in  Appendix  E. 

Two  example  problems  were  worked  using  the  Nelder-Meade  simplex  search 
routine  to  find  the  minimum  performance  index.  When  the  desired  compensator  order 
was  equal  to  the  LQG  compensator  order,  the  numerical  method  and  the  standard  LQG 
solution  method  produced  the  same  results.  Reducing  the  compensator  from  fourth  order 
to  first  order  in  the  second  example  only  increased  the  performance  index  by  0.08  percent 
(page  5-17). 


1  Refers  to  the  page  in  this  thesis  containing  this  conclusion. 
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The  solution  was  overly  sensitive  to  the  initial  guess  for  higher  order 
compensators  (more  than  six  unknown  variables).  Additionally,  the  MATLAB™ 
Lyapunov  routine  should  be  rewritten  so  that  it  does  not  return  an  error  message  and 
interrupt  the  program  when  the  Lyapunov  solution  is  not  unique.  A  large,  arbitrary 
performance  index  should  be  assigned  instead. 

The  numerical  solution  technique  developed  in  this  chapter  was  satisfactory  for 
lower  order  compensators.  The  solution  was  not  sensitive  to  the  initial  guess  and  the 
computing  time  was  not  significant  (page  5-17). 
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6.  The  Sub-Optimal  Pilot  Model 


General 

The  primary  advantage  of  classical  pilot  models  is  that  their  gain,  lead,  lag,  and 
delay  form  is  based  on  experimental  observations  of  human  pilot  behavior.  Choosing 
values  for  these  variables  and  then  using  them  to  predict  a  Cooper-Harper  rating, 
however,  is  difficult.  The  primary  advantage  of  optimal  pilot  models  is  that  they,  like  the 
pilot,  find  the  control  strategy  that  minimizes  a  combination  of  average  tracking  error 
(performance)  and  control  usage  (workload).  This  weighted  combination,  or 
performance  index,  can  be  directly  related  to  a  predicted  Cooper-Harper  rating  based  on 
statistical  fits  to  historical  data.  Because  they  are  based  on  linear  quadratic  Gaussian 
(LQG)  theory,  however,  the  pilot  describing  functions  predicted  by  optimal  pilot  models 
tend  to  be  high  order  and  therefore  inconsistent  with  experimentally  observed  behavior. 
Additionally,  when  the  constraints  of  the  human  pilot  are  implemented  within  the  LQG 
structure,  the  intuitive  nature  of  the  model  is  lost. 

The  sub-optimal  pilot  model  developed  in  this  chapter  uses  the  numerical  LQG 
solution  method  described  in  Chapter  5  to  restrict  the  optimal  pilot  model  solution  to  the 
classical  pilot  model  form.  While  this  solution  may  be  sub-optimal  in  comparison  to  the 
standard  LQG  solution,  it  more  accurately  models  observed  human  pilot  behavior.  The 
sub-optimal  pilot  model  is  not  presented  as  a  solution  to  the  pilot  modeling  problem.  To 
do  so  would  require  die  examination  of  a  broad  handling  qualities  data  base,  requiring 
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years  of  experience  and  model  development.  Rather,  it  is  a  unique  and  promising 
approach,  offered  as  a  prototype  for  future  study. 

This  chapter  is  divided  into  four  sections.  The  first  section,  Model  Development, 
describes  the  sub-optimal  pilot  model.  The  next  section,  Parameter  Analysis,  examines 
the  model's  input  parameters.  In  the  third  section,  Results,  the  Cooper-Harper  ratings  and 
pilot  describing  functions  predicted  by  the  sub-optimal  pilot  model  are  compared  with 
flight  test  results  as  well  as  the  predictions  of  other  pilot  models.  Finally,  the  conclusions 
of  this  chapter  are  summarized. 

Model  Development 

Model  Structure.  The  sub-optimal  pilot  model  is  based  on  the  classical  pilot 
model  structure.  This  can  be  represented  for  the  pitch  axis  as  shown  in  Figure  6-1 


Figure  6-1.  Classical  Pilot  Model  Structure 


1  The  standard  Laplace  variable,  s,  is  used  in  this  chapter  for  convenience.  However,  all  Laplace 
transform  representations  of  the  human  pilot  used  in  this  thesis  are  strictly  valid  only  in  the  frequency 
domain  with  continuous,  random-like  inputs.  They  should  not  be  used  for  system  responses  to 
deterministic  inputs. 
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where 


4  =  Task  Driving  Noise 

rjm  =  Muscular  Noise 

r)o  =  Observation  Noise 

Kp  =  Pilot  Gain 

Tl  =  Pilot  Lead  Time  Constant 

Tj  =  Pilot  Lag  Time  Constant 

x  =  Pilot  Delay 

Tn  =  Muscular  Time  Constant 


0  =  Pitch  Angle 
0C  =  Pitch  Command 
e  =  Task  Error 
F,  =  Stick  Force 
8,  =  Stick  Deflection 
Yc  -  Aircraft  Dynamics 
Yr  =  Task  Dynamics 


This  classical  pilot  model  structure  is  equivalent  to  that  shown  in  Figure  2-5 
(page  2-10),  where  the  pilot  applies  gain  and  washout  to  the  error  and  error  rate  signals. 
This  structure  is  not  identical,  however,  to  that  used  for  the  development  of  the  numerical 
solution  to  the  general  LQG  problem.  The  numerical  method  developed  in  Chapter  5 
must  therefore  be  adapted  slightly. 

For  the  sub-optimal  pilot  model,  the  aircraft  and  task  forcing  function  dynamics 
will  be  assumed  to  be  single-input-single-output  systems.  As  a  result.  Figure  6-1  can  be 
redrawn  for  the  pitch  axis  as  shown  in  Figure  6-2.  Notice  that  the  motor  noise  was 
removed  from  the  model.  Also,  unlike  the  general  LQG  problem  presented  in  Chapter  5, 
the  disturbance  noise  is  now  injected  at  the  controlled  element  output  so  that  the  error 
vector  (0c-0)  can  be  formed. 

Motor  noise  was  not  included  in  the  sub-optimal  pilot  model  for  two  reasons. 
First,  it  was  assumed  that  the  observation  noise  would  adequately  account  for  its  effects. 
Second,  if  motor  noise  is  added  to  the  control  vector,  the  problem  cannot  be  solved 


6-3 


Figure  6-2.  Sub-Optimal  Pilot  Model  Structure 
directly.  The  sub-optimal  pilot  model  weights  control  and  tracking  error  in  the 
following  performance  index. 

oo 

/=  jleTQe  +  uTRu\dt  ( 

o 

The  numerical  solution  assumes  that  the  system  is  driven  by  random  noises  such  that: 

J-E{eTQe  +  uTRu)  ( 

If  noise  is  added  to  the  control  vector,  u,  it  can  be  written  as: 

u  =  [CKxK  +r\m]  (' 


6' 


Substituting  Equation  6-3  into  the  second  term  in  Equation  6-2  results  in  Equation  6-4. 


E{uTRu}=RE{uTu)  =  RE{(Ckxk  +  x\m)2}  (6-4) 

Because  of  the  single-input-single-output  nature  of  the  sub-optimal  pilot  model,  u  is  a 
scalar,  and  several  simplifying  assumptions  were  made  as  shown.  Equation  6-4  expands 
to: 

E{(Ckxk+y\m)2}  =  E{{Ckxkf}  +£{(ru)2}  +E{2CkxK  •  ti„}  (6-5) 

The  first  term  can  be  computed  as  shown  in  the  previous  chapter.  The  third  term  contains 
cross-correlations  that  can  only  be  solved  for  the  single-input-single-output  case.  Finally, 
the  second  term  drives  she  performance  index  to  infinity.  The  reasons  for  this  are  the 
same  as  for  requiring  a  proper  compensator  (page  5-4).  By  definition,  the  expected  value 
of  the  motor  noise  is  as  shown  in  Equation  6-6. 

=/?o-5(t)  (6-6) 

Because  of  the  delta  function,  the  magnitude  of  Equation  6-6  and  the  performance  index 
that  contains  this  term  is  infinite. 

The  equations  of  state  for  the  sub-optimal  pilot  model  system  in  Figure  6-2  are: 


xc  =  AcXc  +  Bcu 

8  =  CcXc 

(6-7) 

u  =  CKxk 

(6-8) 

Xr  =  ArXr  +B 

Qc  =  CrXr 

(6-9) 

xk  =Ak  xK+BK(r\o  +  e) 

ii 

n 

1 

(6-10) 
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Substituting  and  combining  terms  yields  die  following  state-space  representation  driven 
only  by  white  noises. 


Xc 

xr 

— 

Ac  0  BcCk 

0  At  0 

9 

xc 

Xr 

+ 

’  0  0 

Br  0 

m 

' 

_  Xk  . 

-B/cCc  BfcCr  Ax 

.  xx  . 

.  0  BK 

.  T1  . 

where 


e  =  [-Cc  Cr  0  ]• 


xc 

Xr 

xK 


(6-12) 


The  filter  state  space  matrix,  Ar,  must  have  all  eigenvalues  in  the  open  left  half  plane 
since  xr  is  not  controllable  from  u.  This  system  is  in  the  form: 


x=Ax+Ew 


(6-13) 


The  autocorrelation  matrix,  X,  can  be  found  by  solving  the  following  Lyapunov  equation. 

AX+XAT+EQnET=  0  (6-14) 

where  QK  contains  the  noise  intensities,  E;  and  q,  on  the  diagonals.  The  autocorrelation 
matrix,  X,  is  related  to  the  expected  values  of  the  states,  x,  by  the  following  identity. 

X=E{xTx}  (6-15) 

Performance  Index.  The  performance  index  used  in  the  sub-optimal  pilot  model 
weights  error  and  control  such  that: 

J=E{eTQe}+E{uTRu}  (6-16) 
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Because  this  model  is  restricted  to  single-input-single-output  systems,  this  equation  can 
be  rewritten  as: 

J=  Q  E{eTe)  +R  E{uTu)  (6-17) 

where  Q  and  R  are  the  scalar  weightings.  Using  Equations  6-7  through  6-10,  and 
assuming  the  single-input-single-output  case,  this  becomes: 

/=  Q  E{(Crxr  -  Ccxc )2}  +R-  E{xtkCtkCkxk}  (6-18) 

The  second  term  is  identical  to  that  in  Chapter  5  and  can  be  written  for  the 
single-input-single-output  case  as: 

R-ttc^rCK^X^ vm  (6-19) 

«=iy=i 

where  m  is  the  number  of  compensator  states,  n  is  the  number  of  aircraft  and  filter  states, 
and  A' is  the  autocorrelation  matrix  from  Equation  6-15.  The  first  term  in  Equation  6-18 
is  much  more  complicated,  and  can  only  be  solved  for  the  single-input-single-output  case 
where  Crxr  and  Ccxc  are  scalars.  This  term  can  be  expanded  as  in  Equation  6-20. 

QE{{CvXr-Ccxcy}  =  0£{(Cr*r)2}  +QE{{Ccxc)2} 

-2Q-  E{(CcXc )  •  (Crxr)}  (6-20) 

The  first  two  expected  values  are  similar  to  that  in  Equation  6-19.  They  can  therefore  be 
solved  numerically  using  the  following  summations. 

Zs{(Cr*r)2}  =£i  Cr,  ,  •  Crw  ‘X(p+i)^p+j)  (6-21) 

i=iy=i 
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(6-22) 


E{(CcXc)2}  =i  t  Ccu-CciyX,j 

‘-ij-i 

where /is  the  number  of  filter  states,  p  is  the  number  of  aircraft  states,  and  X  is  the 
autocorrelation  matrix  from  Equation  6-15.  The  final  expected  value  in  Equation  6-20  is 
not  as  simple  to  solve.  However,  by  expanding  the  problem  symbolically  and  grouping 
terms,  the  following  summation  can  be  found. 

£{(Crxr)(CcXc)}  =  £  £  Cr, ,  •  CCtJ  •  X^j  (6-23) 

/-u- 1 

In  this  equation,/is  the  number  of  filter  states,  p  is  the  number  of  aircraft  states  and  X  is 
the  autocorrelation  matrix  from  Equation  6-15.  Equations  6-19  and  6-21  through  6-23 
were  used  to  compute  the  performance  index  in  the  sub-optimal  pilot  model.  This 
numerical  algorithm  is  validated  in  Appendix  F. 

Pilot  Model  Form,  Due  to  the  unique  nature  of  the  numerical  approach,  the  LQG 
solution  can  be  restricted  to  a  desired  compensator  form.  Thus,  the  muscular  lag  and 
pilot  delay  can  be  modeled  directly  and  intuitively. 

The  desired  compensator  for  this  model  is: 

c(l)(?  +  c(2))  2-t s 

(s  +  c(3»  (7Vj+1)  2  +  ts  1  ' 

The  vector,  c,  contains  the  three  parameters  the  numerical  routine  will  search  to 
minimize  the  performance  index.  The  first  element  in  the  vector  is  pilot  gain,  the  second 
is  lead,  and  the  is  third  lag.  T„  is  the  muscular  time  constant  and  x  is  the  pilot  delay.  In 
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the  sub-optimal  pilot  model,  the  muscular  lag  time  constant,  TN,  can  be  selected  by  the 
model  user,  but  it  is  not  varied  in  the  numerical  search  algorithm. 

The  pilot  delay  is  modeled  by  a  first  order  Pade  approximation.  Like  the 
muscular  lag  time  constant,  the  delay  is  selected  by  the  model  user,  but  it  is  not  varied  in 
the  numerical  search  algorithm.  A  second  order  Pade  approximation  did  not  affect  the 
results  significantly.  For  the  integrator  example  described  in  Chapter  3  (page  3-8),  the 
performance  index  increased  from  4.6309  to  4.6326  when  a  second  order  Pade 
approximation  was  used  instead  of  a  first.  This  0.037  percent  increase  was  both  typical 
and  insignificant. 

As  a  final  note,  the  numerical  search  algorithm  did  not  consistently  converge 
when  the  pilot  model  was  written  exactly  as  in  Equation  6-24.  Instead,  the  denominator 
had  to  be  expanded  as  shown  in  Equation  6-25. 

- — c(l)s  +  c(2) - 2^j s  (6-25 

S>+(c(3)  +  5L).J+c(3).£  2  +  ” 

Performance  Index  Weightings.  The  standard  optimal  pilot  model  includes 
control  rate  instead  of  control  in  the  performance  index.  It  then  iterates  the  control  rate 
weighting  until  the  desired  muscular  time  constant  is  obtained.  This  establishes  the  ratio 
between  the  weightings  as  shown  in  Chapter  3,  and  makes  the  model's  results 
independent  of  these  weightings  (see  page  3-27  of  this  thesis).  The  sub-optimal  pilot 
model  directly  includes  the  muscular  system  in  the  compensator  form,  leaving  the 
weightings,  Q  and  R,  free  to  be  chosen  by  the  user. 
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Noise  Intensities.  Intensities  must  be  selected  for  both  the  forcing  function  and 
measurement  noises.  Unit  intensity  driving  noise  should  be  used  for  the  task  forcing 
function,  .  The  observation  noise  intensity,  Vn ,  is  iteratively  determined  so  that  the 
noise  ratio,  p,  input  by  the  user  is  obtained. 

The  Cooper-Harper  rating  formula  is  normalized  for  unit  intensity  forcing 
function  noise.  The  observation  noise  intensity,  Vn ,  is  determined  by  solving  the 
numerical  LQG  problem  repeatedly  until  the  desired  noise  ratio  is  achieved.  This  is 
accomplished  using  the  following  relationship. 

Ft,  =  pital  (6-26) 

In  this  equation,  is  the  observation  noise  intensity,  p  is  the  desired  noise  ratio,  and  ce2 
is  the  root  mean  square  magnitude  of  the  task  error  (25:11-14).  In  the  sub-optimal  pilot 
model  structure,  ae2  is  the  first  term  in  Equation  6-18  and  is  found  using  Equations  6-21 
through  6-23.  A  simple  binary  search,  limited  to  five  iterations,  is  used. 

Cooper-Harper  Rating  Prediction.  It  would  be  advantageous  to  use  the  same 
ratings  prediction  formula  used  in  the  STI  optimal  pilot  model,  Equation  3-2  (page  3-8). 
This  formula  was  based  on  the  analysis  of  a  wide  range  of  aircraft  dynamics  and 
normalizes  the  predicted  rating  for  task  intensity.  This  equation  can  not  be  applied  to  the 
sub-optimal  pilot  model  for  two  reasons.  First,  Equation  3-2  was  formulated  for  an 
performance  index  that  weights  control  rate  and  error.  The  sub-optimal  pilot  model 
weights  control  usage.  Second,  the  nature  of  the  two  models  is  quite  different.  The 
sub-optimal  pilot  model  solution  is  restricted  to  a  compensator  form  with  only  three 
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variables  available  for  optimization.  This  produces  a  much  tighter  spread  in  the 
performance  index  values  than  will  occur  in  the  STI  optimal  pilot  model.  Based  on  the 
limited  evaluation  of  the  sub-optimal  pilot  model  conducted  for  this  thesis  the  following 
equations  tend  to  produce  the  best  ratings  correlation.  Note  that  these  equations  are  not 
normalized  for  task  intensity  or  bandwidth.  They  are  based  on  a  regression  analysis  of 
the  dynamics  analyzed  for  this  thesis,  and  are  proposed  as  a  course  measure  to  enhance 
the  parameter  analysis  conducted  later  in  this  chapter. 

Pitch  Rating  =  -30  +  24 1  •  log  10(J) 

Roll  Rating  =  -\3  +  117  ■  log10(J)  ' 

Flow  Diagram.  The  equations  and  computational  flow  of  the  sub-optimal  pilot 
model  are  summarized  in  the  flow  diagram  of  Figure  6-3  on  the  next  page.  The  user 
inputs  the  aircraft  and  task  forcing  functions  along  with  the  error  weighting,  Q,  the 
control  weighting,  R,  the  observation  noise  ratio,  p,  the  muscular  time  constant,  TN,  and 
the  pilot  delay,  td  .  The  sub-optimal  pilot  model  sets  the  task  driving  noise  intensity  at 
one,  the  initial  observation  noise  intensity  at  0.5,  and  all  initial  compensator  coefficients 
to  one.  The  minimum  performance  index  is  then  found  using  the  methods  described 
previously  in  this  chapter.  The  relationship  between  the  observation  noise  intensity  and 
the  desired  noise  ratio  is  examined  using  Equation  6-26,  and  the  process  repeated  (up  to 
five  times)  until  adequate  convergence  is  obtained.  Next,  the  predicted  Cooper-Harper 
rating  is  computed  using  Equation  6-27.  Finally,  the  resulting  pilot  describing  function, 
performance  index,  and  predicted  Cooper-Harper  rating  are  displayed.  The  three 
routines,  or  .to  files,  that  implement  this  model  are  presented  in  Appendix  G.  The  first 
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Where: 

Ac,Bc,Cc  -  Aircraft  Dynamics 
Ag,Bg»Cg  ”  Thsk  Dynamics 
Q  -  Error  Weighting 
R  -  Control  Weighting 
p  “  Observation  Noise  Ratio 
Tn  -  Muscular  Time  Constant 
td»  Pilot  Delay 
K  “  Pilot  Describing  Function 
J  -  Performance  Index 
V  =  Observation  Noise  Intensity 
a2,  =  RMS  Error  Magnitude 


Input: 

Ac,  Be,  Cc,  Ag,  Bg,  Cg 

Q’R’PJn.Td 


Search 


^i  ^min 


Compute  J 


Binary  Search 


Ft,  -  p7tOg 

ZO.ldB 


Output 

Cmin,  J,  Rating 


Figure  6-3.  Sub-Optimal  Pilot  Model  Flow  Diagram 


file,  SOPM.M,  is  the  master  routine.  It  conducts  the  iterations  and  calls  the  other  two 
files  as  necessary.  The  second  file,  PIJSOPM.M,  finds  the  performance  index  value  for 
the  given  compensator,  aircraft,  and  pilot  dynamics.  The  final  file,  RMS_SOPM,  is 
identical  to  PIJSOPM  except  that  in  addition  to  the  performance  index,  it  also  passes 
back  the  root  mean  square  (RMS)  value  of  the  error  for  noise  ratio  computations.  This 
file  is  used  only  after  the  optimal  compensator  is  determined. 


Parameter  Analysis 


The  user  must  choose  several  parameters  when  using  the  sub-optimal  pilot  model. 

1.  Aircraft  Model  Order 

2.  Task  Forcing  Function  (Ag,  Bg,  Cg) 

3.  Error  Weighting  (Q) 

4.  Control  Weighting  (R) 

5.  Observation  Noise  Ratio  (p) 

6.  Muscular  Time  Constant  (TN ) 

7.  Pilot  Delay  (xD) 

The  data  gathered  in  the  following  analysis  indicates  that  the  values  in  Table  6-1  produce 
the  best  results. 


Table  6-1 

Recommended  Parameters  for  the  Sub-Optimal  Pilot  Model 


Parameter 

Symbol 

Recommended  Value 

Aircraft  Model  Order 

Any 

Task  Forcing  Function 

Ag,  Bg,  Cg 

J2 

0.25  s2  +  ^-s+\ 

Error  Weighting 

Q 

1 

Control  Weighting 

R 

4.5  (pitch);  1  (roll) 

Observation  Noise  Ratio 

P 

0.01  (-20  dB) 

Muscular  Time  Constant 

T 

an 

0.115 

Pilot  Delay 

0.37 
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The  reasons  for  each  of  these  choices  are  given  in  the  following  discussion.  Except 
where  noted,  each  of  these  parameters  was  evaluated  individually  with  the  others  held  at 
the  values  shown  in  this  table.  This  parameter  analysis  is  not  meant  to  fine  tune  the 
model  for  use.  The  data  base  evaluated  was  much  to  small  for  this,  and  the  model  is  only 
a  prototype  for  future  study.  Rather,  this  analysis  is  offered  for  the  insights  it  provides. 

The  following  dynamics  and  Cooper-Harper  ratings,  shown  in  Table  6-2,  were 
used  for  this  parameter  analysis.  These  values  are  from  the  in-flight  sum-of-sines 
tracking  task  described  in  Chapter  4. 


Longitudinal  Dynamics:  Lateral  Dynamics: 


Case  1: 

20(5+  1.8)e-°041 
s(s2  +  8.45  +  36) 

Case  A: 

2.5e~°OAs 

5(5 +  2.5) 

(6-28) 

Case  2: 

20(5+  1.8)e~°04i 

5(5 2  +  4.85  +  36) 

Case  B: 

g-0.04s 

5(5+  1) 

(6-29) 

Case  3: 

20(5+  1.8)e~°'24s 
s(s2  +  8.45  +  36) 

Case  C: 

2.5e~°24s 

5(5 +  2.5) 

(6-30) 

Case  4: 

20(5+  1.8)fT°-24s 
s(s2  +4.85  +  36) 

Case  D: 

e-0.24 s 

5(5+  1) 

(6-31) 

Table  6-2 

Cooper-Harper  Ratings  Used  for 
Sub-Optimal  Pilot  Model  Parameter  Analysis 


r 

Case  1 

Case  2 

Case  3 

Case  B 

Case  C 

Case  D 

Rating 

2 

2 

5 

5 

1 

2 

4 

6 
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Aircraft  Model  Order.  Any  aircraft  model  order  can  be  used.  Unlike  the  ST1 
optimal  pilot  model  discussed  in  Chapter  3,  the  order  of  the  pilot  describing  function 
predicted  by  the  sub-optimal  pilot  model  does  not  depend  on  the  aircraft  model  order. 
Additionally,  increasing  the  aircraft  model  order  only  increases  the  size  of  the  Lyapunov 
matrices,  and  does  not  significantly  increase  the  computing  time. 

Some  higher  order  dynamics  are  problematic,  however.  Full  order  dynamics,  by 
their  nature,  have  poorly  scaled  state  space  representations.  For,  example  the  numerator 
for  Case  1  has  a  numerator  gain  of  over  200,000  (refer  to  Equation  C-3),  but  the  first 
term  in  the  denominator  has  unity  gain.  The  numerical  search  routine  did  not  converge 
for  the  full  order  dynamics  evaluated  in  the  flight  test  described  in  Chapter  4.  However, 
the  short  period  approximations  of  Equation  6-28  through  6-31  did  not  present  any 
special  problems. 

Task  Forcing  Function.  When  comparing  model  predictions  with  flight  test 
results,  a  task  forcing  function  consistent  with  the  actual  task  should  be  used.  Otherwise, 
the  following  second  order  Butterworth  filter  with  a  break  frequency,  o>b,  of  two  radians 
per  second  is  recommended. 


Yr(s)  = 


J2 


(6-32) 


This  break  frequency  is  consistent  with  that  normally  used  in  airborne  compensatory 
tracking  tasks.  The  square  root  of  two  in  the  numerator  makes  the  filter  bandwidth  (zero 
dB  crossover)  equal  to  the  filter  break  frequency. 
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As  the  break  frequency  was  moved  from  two,  the  correlation  between  the 
performance  indices  and  the  actual  aircraft  handling  qualities  degraded  rapidly.  For 


example  the  effects  of  decreasing  the  break  frequency  to  0.4  radians  per  second  is  shown 
in  Table  6-3. 


Table  6-3 

The  Effects  of  Task  Break  Frequency  on  Sub-Optimal  Pilot  Model  Results 


■H 

Break 

Frequency 

Performance 

Index 

Minimum  Compensator 
Values1 

Case  1 

2 

1.3637 

0.879,  0,  1.826 

Case  2 

2 

1.3473 

1.218,  0,  1.506 

Case  3 

2 

1.3885 

0.9441,  0,  2.600 

Case  4 

_ 2 

_ 2 

_ 2 

Case  1 

0.4 

0.1620 

3.820,  0,  0.297 

Case  2 

0.4 

0.1544 

2.209,  85.212,  23.546 

Case  3 

0.4 

0.1374 

3.807,0  0.281 

Case  4 

0.4 

0.1333 

3.809,  0,  0.263 

1  Where  the  pilot  is  in  the  form:  (cl*s  +  c2)  exp  (-tau*s)  /  (s  +  c3)*(s  +1/Tn) 

2  A  "matrix  poorly  conditioned"  warning  was  given  while  evaluating  this  case. 


The  performance  indices  correlate  well  with  the  Cooper-Harper  ratings  gathered  during 
flight  test  for  a  break  frequency  of  two.  Notice  that  the  trend  is  reversed  when  the  lower 
break  frequency  is  used.  Such  behavior  makes  this  model  sensitive  to  the  task  and  is 
undesirable. 


Weightings.  It  was  concluded  from  the  flight  test  described  in  Chapter  4,  that 
weighting  stick  deflection  and  task  error  produced  the  best  correlation  with  in-flight 
Cooper-Harper  ratings.  Values  for  the  correlation  coefficients  were  computed  and  are 
presented  in  Table  6-4.  As  shown  in  this  table,  the  control  weighting,  R,  should  be  about 
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4.5  times  the  error  weighting,  Q,  for  pitch  tasks.  The  weightings  should  be  roughly  equal 
for  roll  tasks1. 


Table  6-4 

Sub-Optimal  Pilot  Model  Weighting  Coefficients 


Pitch  Axis 

Roll  Axis 

Tracking  Error  Weighting,  Q 

7.2 

1.7 

Stick  Deflection  Weighting,  R 

32.5 

2.2 

These  values  also  produced  the  best  results  in  the  sub-optimal  pilot  model.  Each 
of  the  pitch  cases  were  evaluated  using  a  stick  deflection  weighting,  R,  of  4.5  and  again 
with  a  weighting  of  one.  The  results  of  this  analysis  are  presented  in  Table  6-5. 

Table  6-5 

The  Effects  of  Control  Weighting  on  Sub-Optimal  Pilot  Model  Results 


HI 

Control 

Weighting 

Performance 

Index 

Minimum  Compensator 
Values1 

Case  1 

4.5 

1.3637 

0.879,  0, 1.826 

Case  2 

4.5 

1.3473 

1.218,  0,  1.506 

Case  3 

4.5 

1.3885 

0.9441,0,  2.600 

Case  4 

—2 

_ 2 

_ 2 

Case  1 

1 

1.2530 

4.064, 0,  2.2035 

Case  2 

1 

1.2154 

3.869, 0,  1.712 

Case  3 

1 

1.2375 

0.950,  0  1.306 

Case  4 

1 

1.3038 

3.860, 0, 2.7536 

1  Where  the  pilot  is  in  the  form:  (cl  *s  +  c2)  exp  (-tau*s)  /  (s  +  c3)*(s  +1/Tn). 

1  A  "matrix  poorly  conditioned"  warning  was  given  while  evaluating  this  case. 


Note  that  the  performance  index  did  not  correlate  as  well  with  the  handling 
qualities  when  the  lower  weighting  was  used.  When  a  weighting  of  one  was  used,  Case  3 


1  The  ratio  between  the  weightings,  Q  and  R,  is  the  salient  issue.  Any  constant  multiplier  of  the  two 
weightings  can  be  factored  to  the  front  of  the  performance  index. 
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had  better  predicted  handling  qualities  (a  lowo*  performance  index)  than  Case  1 .  In  the 
actual  tracking  task,  Case  3  was  given  a  Cooper-Harper  rating  of  5  and  Case  1  was  given 
a  2. 


Observation  Noise  Ratio.  The  observation  noise  ratio  was  experimentally 
measured  at  0.01,  or  -20  dB  (25:12).  Note  that  power  spectra  dB  (10  log10)  is  normally 
used  when  referring  to  this  value.  The  sub-optimal  pilot  model  results  were  evaluated  for 
noise  ratios  of  0.316  (-15  dB)  and  0.0032  (-25  dB).  The  correlation  between  the 
performance  indices  and  the  actual  aircraft  handling  qualities  was  not  affected  by  the 
different  values  of  noise  ratios.  However,  as  the  observation  noise  ratio  increased,  the 
magnitude  of  the  performance  indices  increased.  The  results  of  some  representative  runs 
are  presented  for  the  roll  axis  in  Table  6-6. 

Table  6-6 

The  Effects  of  Observation  Noise  Ratio  on  Sub-Optimal  Pilot  Model  Results 


■m 

Noise 

Ratio 

— 

Performance 

Index 

Minimum  Compensator 
Values1 

Case  A 

-20  dB 

1.3890 

2.179,  0,3.539 

Case  B 

-20  dB 

1.4070 

1.261,0,3.905 

Case  C 

-20  dB 

1.4046 

1.564,  0, 4.677 

CaseD 

-20  dB 

1.4120 

0.826,  0,  5.387 

Case  A 

-15  dB 

1.3987 

1.096,  0, 2.283 

Case  B 

-15  dB 

1.4101 

0.609,  0, 2.592 

Case  C 

-15  dB 

1.4089 

0.728,  0,  3.168 

Case  D 

-15  dB 

1.4131 

0.366, 0,  3.839 

Case  A 

-25  dB 

1.3795 

4.036,  0,  5.715 

Case  B 

-25  dB 

1.4038 

2.411,0,  6.189 

Case  C 

-25  dB 

1.3998 

3.097, 0,  7.345 

CaseD 

-25  dB 

1.4107 

1.695,  0,  8.182 

1  Where  the  pilot  is  in  the  form:  (cl*s  +  c2)  exp  (-tau*s)  /  (s  +  c3)*(s  +1/Tn). 
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Muscular  Time  Constant  The  muscular  time  constant  was  experimentally 
measured  by  McRuer  at  values  between  about  0.08  and  0.12  (15:171).  The  sub-optimal 
pilot  model  results  were  evaluated  using  the  extremes  of  these  values.  The  data  used  for 
this  evaluation  are  presented  in  Table  6-7. 

Table  6-7 

The  Effects  of  Muscular  Time  Constant  on  Sub-Optimal  Pilot  Model  Results 


■u 

Time 

Constant 

Performance 

Index 

Minimum  Compensator 
Values' 

Case  1 

0.115 

1.3637 

0.879,  0,  1.826 

Case  2 

0.115 

1.3473 

1.218,  0,  1.506 

Case  3 

0.115 

1.3885 

0.9441,0, 2.600 

Case  4 

— 

_ 2 

_ 2 

Case  1 

0.083 

1.3641 

0.332, 0, 1.353 

Case  2 

0.083 

1.3472 

1.478,  0,  1.188 

Case  3 

0.083 

1.3889 

1.086,  0,  2.039 

Case  4 

0.083 

1.3790 

1.222,  0,  1.839 

1  Where  the  pilot  is  in  the  form:  (cl*s  +  c2)  exp  (-tau*s)  /  (s  +  c3)*(s  +1/Tn). 

2  A  "matrix  poorly  conditioned"  warning  was  given  while  evaluating  this  case. 


As  shown  in  this  table,  changing  the  muscular  time  constant  had  almost  no  effect 
on  the  performance  index  value  or  the  predicted  pilot  describing  function.  A  value  of 
0.1 15  was  used  for  consistency  with  past  research  (consult  the  discussion  in  Chapter  2  of 
this  thesis  (page  2-8). 

Pilot  Delay.  A  pilot  delay  of  0.35  should  be  used  for  stick  displacement 
command  systems,  and  0.25  seconds  should  be  used  for  force  command  systems.  These 
values  are  based  on  the  flight  test  results  presented  in  Chapter  4  (page  4-13). 
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The  sub-optimal  pilot  model  adds  the  aircraft  delay  to  the  pilot  delay  before 
minimizing  the  performance  index.  Thus,  the  effects  of  altering  pilot  delay  values  can  be 
seen  by  comparing  Cases  1  and  3  or  2  and  4.  In  all  cases  evaluated  for  this  thesis, 
increasing  the  delay  increased  the  performance  index. 


Results 

The  sub-optimal  pilot  model  results  produced  using  the  recommended  values  in 
Table  6-1  are  presented  in  Table  6-8. 


Table  6-8 

Sub-Optimal  Pilot  Model  Results 


Performance 

Index 

Predicted 

Rating 

Actual 

Rating 

Minimum  Compensator 
Values' 

Case  1 

1.3637 

2.5 

2 

0.879, 0, 1.826 

Case  2 

1.3473 

1.2 

2 

1.218,  0,  1.506 

Case  3 

1.3885 

4.4 

5 

0.9441,0, 2.600 

Case  4 

1.45022 

8.9 

5 

5.438,  0,  2.94 12 

Case  A 

1.3890 

0.8 

1 

2.179,  0,3.539 

CaseB 

1.4070 

3.7 

2 

1.261,0,3.905 

Case  C 

1.4046 

3.4 

4 

1.564,  0, 4.677 

Case  D 

1.4120 

4.3 

6 

0.826, 0,  5.387 

1  Where  the  pilot  is  in  the  form:  (cl*s  +  c2)  exp  (-tau*s)  /  (s  +  c3)*(s  +1/Tn). 

2  A  "matrix  poorly  conditioned"  warning  was  given  while  evaluating  this  case. 


As  shown  in  this  table,  the  trend  between  performance  index  values  and  the  actual 
handling  qualities  ratings  was  reasonably  consistent.  The  predicted  Cooper-Harper 
ratings  should  not  be  given  too  much  emphasis.  These  values  were  found  using  a 
regression  formula  optimized  for  these  dynamics  cases.  For  this  reason  they  were 
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omitted  from  the  analysis  in  the  previous  sections.  They  are  provided  here  to  confirm  the 
correlation  between  the  performance  index  and  the  actual  Cooper-Harper  rating. 

Bode  plots  of  the  predicted  pilot  describing  functions  are  presented  in  Figures 
G-l  through  G-4.  Notice  that  in  all  cases,  the  optimal  pilot  describing  function  was  one 
with  pure  lead  in  the  numerator.  This  equates  to  pine  error  rate  feedback  with  washout. 
In  other  words  the  error  gain,  Kp, ,  in  Figure  2-5  (page  2-10)  is  zero.  Such  compensation 
delays  the  magnitude  and  phase  roll-offs  to  the  maximum  extent  possible. 

The  pilot  describing  functions  predicted  by  the  sub-optimal  pilot  model  are  by 
design  more  consistent  with  flight  test  results  than  the  pilot  describing  functions 
predicted  by  the  STI  optimal  pilot  model.  The  STI  optimal  pilot  model  is  much  more 
accurate  predictor  of  Cooper-Harper  ratings,  however.  There  is  a  strong  correlation 
between  the  sub-optimal  pilot  model  performance  index  and  the  actual  Cooper-Harper 
ratings,  but  this  model  lacks  the  maturity  necessary  for  consistent  predictions. 

Conclusions 

The  sub-optimal  pilot  model  developed  in  this  chapter  uses  the  numerical  LQG 
solution  method  described  in  Chapter  5  to  restrict  the  optimal  pilot  model  solution  to  the 
classical  pilot  model  form.  This  model  minimizes  a  performance  index  consisting  of 
error  and  control  usage  and  is  restricted  to  single  axis  use  due  to  the  assumptions 
necessary  to  numerically  compute  the  performance  index  value  (page  6-5)1.  The 
numerical  algorithm  was  validated  in  Appendix  F.  Additionally,  motor  noise  was  not 
modeled  because  it  would  drive  the  value  of  the  performance  index  to  infinity  (page  6-3). 
1  Refers  to  the  page  in  this  thesis  containing  this  conclusion. 
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The  sub-optimal  pilot  model  solution  was  restricted  to  a  gain,  lead,  lag,  delay,  and 
neuromuscular  lag.  The  user  inputs  the  desired  muscular  time  constant  and  delay  value 
(modeled  by  a  first  order  Pade  approximation).  The  model  then  searches  all  values  of 
gain,  lead,  and  lag  to  minimize  the  performance  index.  The  model  is  solved  iteratively 
until  the  desired  observation  noise  ratio  is  achieved.  Finally,  the  predicted  pilot 
describing  function,  Cooper-Harper  rating,  and  minimum  performance  index  are 
displayed. 

A  brief  parameter  analysis  was  performed  in  an  effort  to  gain  insight  into  the 
sub-optimal  pilot  model.  Through  this  analysis,  the  recommended  input  parameters 
presented  in  Table  6-1  (6-13)  were  derived.  More  importantly  this  analysis  lead  to 

several  conclusions  concerning  the  model. 

1 .  Any  aircraft  model  order  can  be  used  as  long  as  the  state  space  representation 
is  not  poorly  scaled.  Neither  the  predicted  pilot  describing  function  nor  the 
required  computing  time  are  affected  by  the  model  order  (page  6-15). 

2.  The  task  forcing  function  should  be  modeled  by  a  second  order  Butterworth 
filter  with  a  break  frequency  of  two.  As  the  break  frequency  was  moved 
from  this  value,  the  correlation  between  the  performance  index  and  the 
predicted  Cooper-Harper  rating  was  diminished  (6-16). 

3.  The  performance  index  weightings  estimated  from  the  flight  test  data 
produced  the  best  results  (6-17). 

4.  The  correlation  between  the  performance  indices  and  the  actual  aircraft 
handling  qualities  was  not  affected  by  changes  in  the  noise  ratios.  However, 
as  the  observation  noise  ratio  increased,  the  magnitude  of  the  performance 
indices  increased  (6-18). 

5.  The  value  of  the  muscular  time  constant  had  no  significant  effect  on  the 
sub-optimal  pilot  model  results  (6-19). 

Finally,  the  results  of  the  sub-optimal  pilot  model  were  compared  to  those  of  STI 
optimal  pilot  model.  By  design,  the  describing  functions  predicted  by  the  sub-optimal 


pilot  model  were  more  consistent  with  flight  test  than  those  of  the  STI  optimal  pilot 
model.  The  STI  optimal  pilot  model  was  a  much  more  accurate  predictor  of 
Cooper-Harper  ratings,  however.  The  sub-optimal  pilot  model  is  a  unique  and  promising 
approach  to  pilot  modeling  and  warrants  further  study. 
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7.  Conclusions  and  Recommendations 

The  most  significant  assumption  made  throughout  this  thesis  was  that  the  pilot 
can  be  modeled  as  a  linear  element.  As  discussed  in  Chapter  2,  this  assumption  is  only 
valid  when  the  task  is  random  appearing  and  within  the  capabilities  of  the  pilot.  Under 
these  conditions,  past  experiments  found  that  the  human  pilot  can  be  modeled  as  a  gain, 
lead,  lag,  delay,  and  a  first  order  muscular  lag  (Reference  15). 

There  are  currently  three  broad  categories  of  pilot  models.  The  first  category, 
open  loop  models ,  use  the  aircraft  response  to  open  loop  commands  to  predict  handling 
qualities.  Because  these  models  make  no  attempt  to  directly  model  the  human  pilot,  they 
are  relatively  simple  to  develop  and  are  the  most  commonly  used.  The  second  category, 
classical  pilot  models,  model  the  pilot  as  a  gain,  lead,  lag,  and  delay.  This  form  is  based 
on  experimental  observations,  but  using  these  models  to  predict  Cooper-Harper  ratings 
remains  a  difficult  task.  The  final  type  of  model,  optimal  pilot  models,  model  the  pilot  as 
an  optimal  regulator  and  estimator.  Relating  the  predicted  pilot  describing  function  to  a 
predicted  Cooper-Harper  rating  is  more  straight  forward  with  these  model.  Their 
implementation  is  difficult,  however,  and  the  pilot  describing  functions  they  predict  are 
generally  high  order,  and  therefore  not  consistent  with  observed  human  behavior. 

An  optimal  pilot  model  developed  by  Systems  Technology,  Incorporated,  (STI) 
was  analyzed  in  detail.  This  model  was  chosen  for  three  reasons.  First,  it  incorporates 
nearly  every  important  aspect  of  other  optimal  pilot  models,  making  it  a  good  candidate 
for  study.  Second  it  can  be  implemented  on  the  personal  computer  and  is  therefore 
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widely  available.  Finally,  this  model  has  had  some  success  in  predicting  Cooper-Harper 
ratings,  but  lacks  parameter  selection  guidance. 

The  analysis  of  this  model  included  a  model  overview,  an  integrator  example,  and 
a  detailed  parameter  analysis.  The  integrator  example  was  worked,  step-by-step,  to 
clarify  the  model's  logic.  This  example  paralleled  an  example  presented  in  the  STI 
documentation  (Reference  25).  The  handling  qualities  rating  predicted  for  this  example 
was  not  realistic  (negative),  illustrating  the  need  for  proper  selection  of  the  model's 
parameters. 

The  user  must  select  nearly  twenty  parameters  when  running  the  STI  optimal 
pilot  model.  Most  of  the  selections  are  straight  forward,  and  some  have  no  significant 
effect  on  the  model  results.  Several  important  recommendations  were  made,  however. 
First,  the  lowest  feasible  aircraft  model  order  should  be  used.  The  STI  model  predicts  a 
pilot  describing  function  of  order  2n+5  where  n  is  the  number  of  aircraft  and  filter  states. 
Second,  the  task  bandwidth  should  be  selected  by  running  the  model  for  a  range  of 
bandwidths  and  selecting  the  one  that  produces  the  worst  Cooper-Harper  rating.  If  the 
generally  accepted  flight  test  bandwidth  of  2  radians  per  second  is  used,  the  model  does 
not  produce  satisfactory  results.  Third,  the  filter  should  be  augmented  to  the  plant 
output.  Finally,  the  experimentally  estimated  observation  noise  ratios  of  -20  decibels 
should  be  used.  For  the  best  results,  however,  the  motor  noise  ratio  should  be  increased 
from  its  experimentally  based  value  of  -25  decibels  to  -20  decibels. 

When  the  proper  model  parameters  were  used,  the  STI  optimal  pilot  model  was 
moderately  successful  at  predicting  Cooper-Harper  ratings,  both  for  the  LAMARS  cases 
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analyzed  in  Chapter  3  and  the  flight  dynamics  cases  described  in  Chapter  4. 
Unfortunately,  the  STI  model  predicted  a  fifteenth  order  pilot  when  short  period 
dynamics  were  analyzed.  Bode  plots  of  these  pilot  describing  functions  had  deep 
notches,  indicative  of  this  high  order  compensation.  Due  to  their  high  order,  the  pilot 
describing  functions  predicted  by  the  STI  model  were  not  consistent  with  the  classical 
pilot  model  form  or  observed  human  pilot  behavior. 

A  limited  evaluation  of  human  pilot  response  was  sponsored  in  support  of  this 
thesis  by  the  Air  Force  Flight  Dynamics  Directorate.  Five  sorties  were  flown  in  the 
Calspan  variable  stability  Lear  II  aircraft.  Ground  simulations  in  Lear  II  were  also 
performed.  Four  different  pitch  and  four  different  roll  axis  dynamics  were  evaluated 
using  three  different  tracking  tasks.  Primary  pilot  response  parameters  were  recorded 
and  examined  using  statistical  and  Fourier  transform  analysis  in  an  attempt  to  provide 
insight  into  human  pilot  behavior. 

The  dynamics  evaluated  during  the  flight  test  represented  a  broad  range  of 
handling  qualities  as  evidenced  by  the  assigned  Cooper-Harper  ratings.  Additionally,  the 
variability  in  these  ratings  was  acceptably  low. 

An  analysis  of  the  discrete  tracking  task  time  histories  revealed  that  the  pilot 
delay  between  task  command  and  stick  force  (0.27  seconds)  was  consistent  with  previous 
experimental  estimations.  The  stick  deflection  lagged  stick  force  by  0.1  seconds  in  all 
cases,  however,  due  to  the  lag  effects  of  the  stick  dynamics.  This  produced  a  total  delay 
between  task  command  and  stick  deflection  of  0.37  seconds,  well  above  that  normally 
used  in  pilot  model  analysis. 
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A  regression  analysis  of  the  flight  test  data  was  conducted  in  an  attempt  to 
evaluate  the  validity  of  two  optimal  pilot  model  weighting  schemes.  First,  a  regression 
analysis  of  root  mean  square  (RMS)  tracking  error  and  normalized  RMS  stick  deflection 
to  Cooper-Harper  rating  was  performed.  Second  a  regression  analysis  of  RMS  tracking 
error  and  RMS  stick  deflection  rate  to  Cooper-Harper  rating  was  performed.  The  task 
error  and  normalized  stick  deflection  weighting  scheme  produced  a  much  higher  rating 
correlation  for  the  airborne  data. 

A  frequency  response  analysis  of  task  error  to  stick  deflection  was  conducted  in 
an  attempt  to  provide  insight  into  human  pilot  response.  This  analysis  revealed  that  the 
pilot?  -v  i  not  exhibit  higher  order  behavior.  The  frequency  responses  were  consistent 
with  the  classical  gain,  lead,  and  lag  form  except  for  the  presence  of  large  amounts  of 
pure  phase  lead  at  higher  frequencies.  In  all  cases  the  pilot  acted  so  that  the  combined 
pilot-aircraft  system  resembled  and  integrator  in  the  crossover  region.  When  the  aircraft 
command  path  gain  was  doubled,  the  pilot  reduced  his  gain  so  that  the  response  of  the 
combined  pilot-aircraft  system  was  not  significantly  changed. 

A  numerical  solution  to  the  linear  quadratic  Gaussian  (LQG)  problem  was 
derived  in  Chapter  5.  This  solution  allows  the  compensator  form  to  be  predetermined. 
This  method  may  be  applicable  for  any  situation  when  reduced  order  compensation  is 
desired. 

The  numerical  method  assumes  the  aircraft,  filter,  and  compensator  dynamics  are 
proper.  Further,  the  LQG  weighting  matrices  must  be  diagonal  and  the  system  must  only 
be  driven  by  random  noises.  The  value  of  the  standard  LQG  performance  index  is 
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determined  numerically,  using  the  covariance  matrix  and  a  series  of  summations  for  any 
compensator,  controlled  element,  and  filter.  A  Nelder-Meade  simplex  search  is  then  used 
to  find  the  compensator  coefficients  that  minimize  the  performance  index. 

Two  examples  were  worked  using  this  numerical  approach.  The  first 
demonstrated  that  when  the  predetermined  compensator  order  was  the  same  as  that  of  the 
standard  LQG  solution,  the  two  methods  produced  identical  results.  The  second  example 
demonstrated  the  potential  of  this  method  for  use  in  finding  reduced-order,  sub-optimal 
compensators. 

These  examples  also  revealed  a  deficiency  in  the  numerical  search  routine.  The 
Nelder-Meade  search  routine  was  found  to  be  unsatisfactory  for  complex  compensator 
forms.  When  numerous  parameters  had  to  be  searched  the  success  of  the  routine 
depended  greatly  on  the  initial  guess.  The  method  was  found  to  be  satisfactory  for 
compensators  with  six  or  less  parameters  to  search.  Further,  it  was  recommended  that  the 
MATL  AB™  Lyapunov  solver  be  rewritten  so  that  it  returns  an  arbitrary  and  large 
performance  index  when  the  Lyapunov  solution  is  not  unique,  rather  than  an  error 
message. 

Finally,  the  numerical  LQG  solution  was  used  in  a  sub-optimal  pilot  model 
developed  in  Chapter  6.  This  model  restricted  the  optimal  pilot  model  solution  to  the 
classical  pilot  model  form.  It  was  sub-optimal  in  terms  of  the  standard  LQG 
performance  index  due  to  its  low  order  form,  but  it  was  by  nature  more  consistent  with 
human  pilot  behavior.  This  model  minimized  a  performance  index  consisting  of  task 
error  and  control  usage. 
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Additional  numerical  summations  were  developed  to  implement  the  output 
disturbance  form  and  error  weighting  of  a  compensatory  tracking  task.  The  model  was 
restricted  to  single  axis  use  due  to  the  assumptions  necessary  to  numerically  compute  the 
performance  index  value.  These  new  summations  were  verified  by  example  in 
Appendix  F. 

The  sub-optimal  pilot  model  solution  was  restricted  to  a  gain,  lead,  lag,  delay,  and 
muscular  lag.  The  user  input  the  desired  muscular  time  constant  and  delay  value.  The 
model  then  searched  all  values  of  gain,  lead,  and  lag  to  minimize  the  performance  index. 
The  model  was  solved  iteratively  until  the  desired  observation  noise  ratio  was  achieved. 
Finally,  the  predicted  pilot  describing  function,  Cooper-Harper  rating,  and  minimum 
performance  index  were  displayed. 

A  brief  parameter  analysis  was  performed  in  an  effort  to  gain  insight  into  the 
sub-optimal  pilot  model.  Through  this  analysis,  recommended  input  parameters  were 
derived.  More  importantly,  this  analysis  lead  to  several  conclusions  concerning  the 
model.  First,  any  aircraft  model  order  could  be  used  so  long  as  the  state  space 
representations  were  properly  scaled.  Second,  a  task  forcing  function  consistent  with  that 
used  in  flight  test  could  be  used.  Third,  the  performance  index  weightings  obtained  from 
the  flight  test  regression  analysis  produced  good  results.  Finally,  in  all  cases  the 
predicted  pilot  describing  function  had  a  free  s  in  the  numerator,  generating  as  much  lead 
at  as  low  frequency  as  possible. 

By  design,  the  describing  functions  predicted  by  the  sub-optimal  pilot  model  were 
more  consistent  with  flight  test  than  those  of  the  STI  optimal  pilot  model.  The  STI 
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optimal  pilot  model  was  a  much  more  accurate  predictor  of  Cooper-Harper  ratings, 
however.  The  sub-optimal  pilot  model  is  a  unique  and  promising  approach  to  pilot  model 
and  warrants  further  analysis. 
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Appendix  A.  Nonstandard  Performance  Indices 


The  optimal  pilot  model  discussed  in  Chapter  3  uses  the  following  nonstandard 
performance  index. 

J- 1  [xT(t)  Qc  x(t)  +  u  T(t)  Rc  u(0]  dt  (A- 1 ) 

o 

Notice  the  weighting  penalty  is  on  control  rate  instead  of  control  usage.  Minimizing  this 
performance  index  is  equivalent  to  augmenting  all  channels  of  the  plant  with  integrators 
and  minimizing  the  standard  performance  index  in  Equation  A-2. 

J=\  [*r(0  Qc  x(i)  +  u  T(t)  Rc  u(0]  dt  (A-2) 

o 

Consider  the  standard  LQG  diagram  in  Figure  A-l  below  with  r,  the  commanded 
input,  equal  to  zero  and  the  weighting  matrices,  Q  and  R,  as  shown. 


Figure  A- 1 .  Standard  Control  Weighting 


Minimizing  the  standard  quadratic  performance  index,  Equation  A-2,  is  equivalent  to 
minimizing  the  two-norm  of  the  output  of  the  two  vectors,  Ru(t)  and  Qx(t),  over  an 
infinite  time  horizon.  If  it  is  desired  to  minimize  RU(i)  instead  of  Ru(t),  the  diagram  can 
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be  modified  as  shown  in  Figure  A-2. 


Figure  A-2.  Control  Rate  Weighting 

The  new  plant,  G*,  is  equivalent  to  the  original  plant  augmented  with 
integrators.  The  control  rate  weighting  matrix,  R,  can  now  be  used  to  establish 
compensator  lags  on  the  different  control  channels. 
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Appendix  B.  Predicted  Describing  Functions 


General 


This  appendix  contains  Bode  plots  of  the  pilot  describing  functions  predicted  by 
the  STI  optimal  pilot  model  during  the  sensitivity  analysis  in  Chapter  3.  Each  figure  also 
presents  the  lower  order  equivalent  system  match  to  the  classical  pilot  model  form. 

The  dynamics  analyzed  were: 


Case  1:  jr- 
o« 


100 

5(5  + 100) 


(B-l) 


Case  2:  |-  =  2^— 125) 
s(s2  +  8^  +  25) 

Case  3:  |- =■  2Q(5  +  L25)  •  e-2* 
s(s2  +  8.9  +  25) 

Case  4:  |_=  20(^ +1.25)  g_.033, 

8e  5(52  + 1.85  +  25) 


Case  5: 


0_  =  20(5+ 1.25)  e_.2j 

s(s2  + 1.85  +  25) 


(B-2) 

(B-3) 

(B-4) 

(B-5) 


These  dynamics  have  the  following  characteristics: 


Table  B-l 

STI  Optimal  Pilot  Model  Evaluation  Dynamics 


Case  1 

Short  Period 
Damping  Ratio, 

Short  Period  Natural 
Frequency,  coro 

Delay,  x 

1 

100 

0 

Case  2 

0.8 

5 

0.033 

Case  3 

0.8 

5 

0.2 

Case  4 

0.18 

5 

0.033 

Case  5 

0.18 

5 

0.2 
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The  STI  optimal  pilot  model  predicted  the  following  pilot  describing  functions': 


_  174(0X.0448)[.707,.4](1.65X5.79X12.63)(12.91)[-.866>  17.31(100)2 
p  1  (0)[.707,  .4]2(6.02X12.57X12.63)[.274, 18.96](45.4X94X100) 


92.6(0)(.032)[.707,  .4](1.4S)[.814,4.56][.8,  5j(5.08XI2.5)2[-.866>  14.9] 
(0)[.707,.4]2(I.25X4.96)[.8,5](12.43X12.46)[.  146, 14.2J[.84, 27.3]  C  ’ 


92.7(0)(.035)[.707,.4](U9)(4.1)[.79  6,4.78][.8,5][-.866,8.66](12.5)2 
(0)[.707,.4]2(1.25)(4.02)[.8,5][.043,11.5](12.36)(12.46)[.81,23.5^  (  ‘  ] 

46.2(0)(.028)[.707)  .4](1.18)[.06,4.4](4.56)[.18, 5](12.5)(12.6)[-.866, 14.9] 
(0)[.707,.4]2(1.25X4.6)[.18,5](12.49X12.52)[.157,12.8][.846,23.3]  1  '  * 


40.4(0)(.028)[.707>  .4](1.04X3.74)[.02S,4.96],  [.18, 5],  [-.866, 8.66](12.S)2 
(0)[.707,.4]2(1.25)(3.77)[.182,5][.082,10.7](12.4)(12.5)[.818, 19.12]  (  '  ’ 


The  following  lower  order  equivalent  system  matches  to  the  classical  pilot  model  form 
were  found: 


i  pi 


lP  2 


h  3 


*  p  4 


rpS 


3.756(5  +  2.921)  __ 

(5+1.546) 

(B- 11) 

1.29(f  +  3.248) 

(5  +  0.973) 

(B-12) 

1.294(1  +  2.253) 

(5  + 1.036) 

(B-13) 

0.347(5  +  5.42)  om, 

(5  +  0.682) 

(B-14) 

0.295(5  +  5.31) 

(5  +  0.765) 

(B-15) 

1  Where  the  brackets  denote  [£,  ©,]  as  in  s2+2C©,+<o11s,  and  the  parentheses  denote  (r)  as  in  s+r. 
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Frequency  (radians  per  second) 


Figure  B-l.  Bode  Magnitude  Plot  of  Case  1 


Predicted  Pilot  Describing  Function 
Low  Order  Equivalent  System  Match 


10° 

Frequency  (radians  per  second) 

Figure  B-2.  Bode  Phase  Plot  of  Case  1 


Figure  B-3.  Bode  Magnitude  Plot  of  Case  2 


Figure  B-4.  Bode  Phase  Plot  of  Case  2 


B 


Frequency  (radians  per  second) 


Figure  B-6.  Bode  Phase  Plot  of  Case  3 
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Predicted  Pilot  Describing  Function 
Low  Older  Equivalent  System  Mate 


10° 

Frequency  (radians  per  second) 


Figure  B-9.  Bode  Magnitude  Plot  of  Case  5 


Frequency  (radians  per  second) 


Figure  B- 10.  Bode  Phase  Plot  of  Case  5 


Program  CC  Macro 


The  following  macro  implements  the  STI  optimal  pilot  model  on  Program  CC  as 
described  in  Chapter  3.  These  commands  use  the  baseline  parameters  listed  in  Table  3-3 
(page  3-21),  where  ycl  is  the  aircraft  transfer  function  and  yw  is  the  task  forcing 
function. 


File  Name:  *.MAC  in  OCM  subdirectory 

@ocmyc  1  ,yc  1  ,yw,p40 

@ocmlqr,p40,l,.l,.08,.001 

@ocmsetup,  1  ,.2, -20, -20, -20, 1,1,1 ,0,0, 1 

@ocmkbf,.l 

@ocmpilot,.2,2 

state 

gep,p24,yp 

cc 

yp=-yp 

yp 
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Appendix  C.  Flight  Test  Information 


General 

This  appendix  presents  supplemental  information  for  Chapter  4,  Flight  Test. 


Dynamics  Implementation 

The  dynamics  described  on  page  4-8  of  this  thesis  were  implemented  as  a  position 
command  system  as  shown  in  Figure  C-l  and  C-2. 


Commanded 
Stick  Deflection 


Figure  C-l.  Feel  System 


rVtmmaivfoH  Aircraft 


Figure  C-2.  Flight  Control  System 


The  elevator  and  aileron  actuator  dynamics  were: 


702 

s2  +  2(.7)(70)  +  702 


(C-l) 


The  longitudinal  and  lateral  stick  dynamics  were: 


162 

j2  +2(.7)(16)  + 162 


C-l 


The  feel  system  characteristics  were: 


Elevator: 

Aileron: 

Stick  Force  Gradient 

6  Ib/in 

4  lb/in 

Stick  Breakout  Force 

0.75  lb 

0.75  lb 

Stick  Force  per  g 

7  lb/g 

Control  Gearing 

8  deg/in 

12  deg/in 

The  following  equations  represent  the  linear  implementation  of  these  dynamics1 
(31:49-50). 


8  _  16*  1  8  _  70*  5.5(1. g) 

F„  [.7,16]  6  [.7,70]  (0)[GP,6] 


8  -  6  70*  5.5(1, 8) 

*«  [-7. 70]  (0)K„,6] 


1 - 16*  1  ,,  ,  702  3.3 

F-  [  7, 16]  4  [.7,70]  mTR) 


jL-12.-72i_._iJ_.,-*, 
*«  [-7,70]  mn) 


(C-3) 

(C-4) 

(C-5) 

(C-6) 


Test  Point  Evaluation  Cards 

The  test  point  evaluation  card,  presented  in  Figure  C-3,  and  the  pilot  induced 
oscillation  (PIO)  rating  scale,  presented  in  Figure  C-4,  were  used  during  the  flight  test 
described  in  AFFTC-TLR-93-41.  They  are  presented  here  for  the  reader's  convenience. 
1  Where  [  7,16]  denotes  [£,  coj  as  in  sJ+2^con+con2  and  (1.8)  denotes  (x)  as  in  (s+x). 
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CASE 

Itera- 

Sortie 

# 

tion# 

# 

HAVE  PILOT  TEST  CARD 


Filename 


PRE-BRIEF 

All  test  points  start  at  15,000  feet  MSL,  250  KIAS 

The  pilot  will  perform  each  task  using  any  pilot  compensation  necessary  to  minimize  the 
average  tracking  task  error. 

DESIRED:  PIO  tendencies  do  not  compromise  tracking  task.  Commanded  attitude 
maintained  within  0.5  degrees  in  pitch  and  5  degrees  in  bank  (measured  at  end  of 
command  bar)  for  50%  of  the  task  except  immediately  following  step  command  changes. 

ADEQUATE:  Commanded  attitude  maintained  within  1  degree  in  pitch  and  10  degrees  in 
bank  (measured  at  end  of  command  bar)  for  50%  of  the  task  except  immediately  following  step 
command  change. 


POST-BRIEF  C-H  PIO 

Rating  Rating 

1.  Assign  PIO  rating  - 

2.  Assign  Cooper-Harper  rating 

3.  Aircraft  response  to  input  (pitch/roll) 

Initial  -  Quick,  Slow,  Sluggish,  etc 
Final  -  Predictable,  Crisp,  etc 

4.  Does  the  level  of  aggressiveness  affect  task  performance  (precision,  accuracy,  etc)? 

5.  Any  special  piloting  techniques/compensation  required? 

6.  Any  undesirable  aircraft  motions  (turbulence,  disorienting)? 

7.  Provide  actual  percentage  performance  to  pilot  Turb  %  in  %  in  C-H 

Rating  Desired  Adequate  Re-Rating 

8.  Review  Cooper-Harper  rating  ^ 


Figure  C-3.  Test  Point  Comment  Card 


Figure  C-4.  Calspan  Pilot  Induced  Oscillation  Rating  Scale 


Software  Validation  Test  Case 

Figure  C-5  presents  the  magnitude  (top  plot),  phase  (middle  plot),  and  coherence 
(bottom  plot)  computed  by  the  frequency  response  analysis  software  as  described  in 
Chapter  4  of  this  thesis  (page  4-10).  Each  of  these  is  plotted  as  a  function  of  frequency. 
Individual  data  points  are  represented  by  asterisks.  The  dotted  lines  represent  95  percent 
confidence  interval  bounds  (7:55). 

Figure  C-6  presents  the  Bode  plot  found  using  the  frequency  response  analysis 
software  along  with  the  actual  Bode  plot  of  the  linear  system.  Except  for  the  lowest 
frequency  phase  point,  the  two  plots  are  nearly  identical. 
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Figure  C-5.  Frequency  Response  Analysis  --  Validation  Case 


Phase  (degrees)  Magnitude,  dB  (decibels 


Frequency,  to  (radians  per  second) 


Frequency,  co  (radians  per  second) 

Figure  C-6.  Bode  Comparison  Plot  —  Validation  Case 

Cooper-Harper  Ratings  and  Pilot  Comments 

The  Cooper-Harper  ratings  and  pilot  comments  for  the  ground  and  airborne 
evaluations  are  summarized  in  Tables  C-l  and  C-2.  Only  the  single  axis,  sum-of-sines 
evaluations  from  AFFTC-TLR-93-41  were  analyzed  in  this  thesis  (7:73-107). 
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Table  C-l 

Cooper-Harper  Ratings  and  Pilot  Comments  —  Ground  Simulation 


Case 

Iteration 

(sortie) 

PIO 

Rating 

C-H 

Rating 

Percent 

Desired 

Percent 

Adequate 

Comments1 

1 

1(1) 

1 

2 

49 

90 

Good  response  in  the  pich  axis.  I  didn't 
have  to  be  overly  aggressive 

1 

2(3) 

1 

3 

52 

91 

A  little  slow  in  pitch.  The  more 
aggresive  I  got  the  better  I  could  track. 

2 

1(2) 

2 

H 

54 

86 

Satisfactory.  Initial  respone  was  quick. 
Tendecny  to  overshoot. 

2 

2(3) 

2 

3 

39 

76 

Easy  to  acquire.  Couple  of  overshoots. 
Slightly  sensitive.  Maybe  too  quick. 

3 

1(1) 

3 

5 

43 

82 

Large  inputs  required.  Slight  PIO.  A 
little  sensitive.  Lead  required. 

3 

2(6) 

2 

H 

39 

81 

Slow  initial  response.  Sluggish. 

Moderate  compensation  required. 

B 

1(1) 

H 

6 

44 

72 

Large  inputs  required.  PIO  Tendency. 

Have  to  come  out  of  the  loop. 

B 

2(5) 

2 

4(5) 

35 

68 

Undesireable  motions.  Slow.  Tend  to 
overshoot.  Unpredictable.  Required  lead. 

B 

1(2) 

1 

1 

80 

99 

Quick.  Predictable.  Fairly  aggressive. 
Pilot  compensation  not  a  factor. 

B 

2(3) 

2 

2 

88 

99 

Slow  response.  Easy  to  reacquire.  No 
undesirable  motions.  Comp,  nota  factor. 

B 

1(4) 

1 

2 

82 

100 

Initial  input  was  good.  No  tendency  to 
overshoot.  Work  load  not  too  high. 

B 

2(5) 

1 

2 

75 

98 

Quick.  Predictable.  Could  be  fairly 
aggressive.  Stayed  in  the  loop. 

C 

1(5) 

2 

H 

63 

91 

Need  to  back  out  a  little.  Not  slow,  but 
not  predictable.  Slight  ocillatory  motion. 

C 

2(6) 

2 

3 

77 

99 

Overshoots.  A  little  compensation 
required.  Tend  to  back  out  of  loop.  Slow. 

D 

1(2) 

3 

5 

61 

87 

Sluggish.  Tendency  to  overshoot.  Large 
lead  input  required.  Backed  out  of  loop. 

D 

2(4) 

2 

3 

77 

96 

Tendency  to  overshoot.  Work  load 
tolerable.  Little  oscillation,  but  not  obj. 

1  These  comments  were  summarized  from  AFFTC-TLR-93-41  (7:83-90). 
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Table  C-2 

Cooper-Harper  Ratings  and  Pilot  Comments  —  Airborne  Evaluation 


Case 

Iteration 

(sortie) 

PIO 

Rating 

C-H 

Rating 

Percent 

Desired 

Percent 

Adequate 

Comments* 

1 

1(1) 

1 

2 

67 

98 

Predictable.  No  PIO  tendency.  Good 
initial  response.  No  real  comp.  used. 

1 

2(2) 

2 

3 

68 

97 

Backed  out  of  the  loop  to  prevent 
oscillation  tendencies. 

1 

3(5) 

1 

2 

69 

96 

Good  initial  response.  Slight  bobble  but 
negligible  deficiencies. 

2 

KD 

1 

1 

53 

93 

Didn't  push  as  hard  as  I  did  on  the 
ground. 

2 

2(2) 

1 

2 

63 

98 

No  PIO  tendency.  Compensation  not  a 
factor.  Backed  out  a  little  bit. 

2 

3(5) 

1 

2 

70 

97 

Good  initial  response.  A  little  bobble. 
Don't  have  to  be  aggressive  at  all. 

3 

1(1) 

H 

5 

45 

86 

Oscillatory  tendency.  Need  to  come  out 
of  the  loop.  PIO  tendency.  Unpredictable. 

3 

2(2) 

H 

5 

48 

97 

Can  track  fairly  aggressively.  Moderately 
objectionable.  Slow  initial  response. 

3 

3(5) 

2 

3 

62 

98 

No  tendency  to  oscillate.  Lagging  with 
my  inputs.  Initial  response  a  bit  slow. 

H 

1(1) 

B 

■ 

47 

89 

PIO  tendency.  Can't  be  too  aggressive. 
Light  on  the  controls.  High  worload. 

H 

2(2) 

2 

B 

48 

80 

Oscillatory  especilally  when  aggressive. 
Less  than  predictable.  Tend  to  overshoot. 

B 

3(5) 

3 

5 

41 

76 

Tend  to  overshoot.  Tend  to  back  out  a 
little  bit.  Osillalions.  Not  predictable. 

1  These  comments  were  summarized  from  AFFTC-TLR-93-41  (7:83-90). 
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Table  C-2  (Continued) 

Cooper-Harper  Ratings  and  Pilot  Comments  —  Airborne  Evaluation 


Case 

Iteration 

(sortie) 

PIO 

Rating 

C-H 

Rating 

Percent 

Desired 

Percent 

Adequate 

Comments1 

D 

KD 

1 

1 

91 

100 

No  tendency  to  overshoot.  Predictable 
response.  Like  the  way  it  handled. 

D 

2(2) 

1 

3 

86 

99 

Fine  tracking  is  simple.  No  oscillation  at 
all.  Response  predictable. 

H 

3(5) 

1 

1 

95 

99 

Can  make  prety  quick  inputs.  Quick. 
Predictable.  No  compensation,  required. 

B 

1(1) 

1 

2 

86 

99 

Initial  response  OK.  Put  in  big  input  then 
take  it  out.  Anticipation  required. 

B 

2(2) 

2 

B 

82 

99 

Fine  tracking  not  a  problem.  Slight 
oscillatory  tendency.  Lacks  predictability. 

B 

3(5) 

1 

2 

90 

99 

Response  initially  good.  Easy  to  track 
target.  Nice  handling  aircraft. 

C 

KD 

2 

3 

74 

100 

If  too  aggressive,  get  more  oscillations. 
Response  is  sluggish  and  unpredictable. 

C 

2(2) 

2 

B 

79 

99 

Can't  be  real  aggressive  in  fine  tracking. 
Slight  delay. 

C 

3(5) 

2 

B 

79 

99 

Fine  tracking  not  too  bad.  Tendency  to 
over-control  with  large  inputs. 

D 

1(1) 

H 

6 

62 

93 

Had  to  back  out  of  the  loop.  Required 
lots  of  lead. 

D 

2(2) 

H 

B 

71 

98 

Practically  flying  this  open  loop.  Initial 
response  was  slow.  Overshoot  tendency. 

D 

3(5) 

3 

B 

77 

98 

Large  lead  input  required.  Oscillations. 
Came  out  of  the  loop.  High  work  load. 

1  These  comments  were  summarized  from  AFFTC-TLR-93-41  (7:83-90). 
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Figure  C-7.  Cooper-Harper  Ratings 


Statistical  Analysis 


The  statistics  in  Tables  C-3  and  C-4  were  computed  using  the  flight  test  data  base 
produced  by  AFFTC-TLR-93-41  and  are  summarized  in  Table  4-2.  Only  the  single  axis 
sum-of-sines  data  were  analyzed  for  this  thesis. 


Table  C-3 

Optimal  Pilot  Model  Statistical  Analysis  Parameters 
—  Ground  Simulation  Single  Axis  Sum-of-Sines  Data 


Case 

Sortie# 

C-H  Rating 

RMS  Error 

NRMS  Stick 

RMS  Stick  Rate 

1 

3 

3 

0.5894 

0.3306 

0.6290 

2 

2 

4 

0.6738 

0.3654 

0.6853 

2 

3 

3 

0.7911 

0.3374 

0.5088 

3 

1 

5 

0.7506 

0.3478 

1.1005 

3 

6 

4 

0.7242 

0.3902 

0.8480 

4 

1 

6 

0.8765 

0.3910 

1.2988 

4 

5 

5 

0.9031 

0.4135 

1.0112 

A 

2 

1 

3.9868 

0.3490 

1.1969 

A 

3 

2 

3.1897 

0.3752 

1.1342 

A 

7 

2 

3.4809 

0.3727 

1.1294 

B 

4 

2 

3.5883 

0.2947 

1.0995 

B 

5 

2 

4.2318 

0.3069 

1.0815 

C 

5 

4 

5.4722 

0.4556 

1.7275 

C 

6 

3 

3.7028 

0.4496 

1.3451 

D 

2 

5 

6.2632 

0.4538 

1.4944 

D 

4 

3 

4.7876 

0.3127 

1.1310 

where: 


RMS(x)  = 


V  n 


(C-7) 


and  n  is  the  number  of  samples  in  the  vector,  x. 


Table  C-4 

Optimal  Pilot  Model  Statistical  Analysis  Parameters 
—  Airborne  Single  Axis  Sum-of-Sines  Data 


Case 

Sortie  # 

C-H  Rating 

RMS  Error 

NRMS  Stick 

RMS  Stick  Rate 

1 

1 

2 

0.4649 

0.3590 

0.6647 

1 

2 

3 

0.4633 

0.3931 

0.6542 

1 

5 

2 

0.4694 

0.3693 

0.7803 

2 

1 

1 

0.6112 

0.3097 

0.7035 

2 

2 

2 

0.4996 

0.3364 

0.7024 

2 

5 

2 

0.4462 

0.3546 

0.6393 

3 

1 

5 

0.6715 

0.3770 

1.0363 

3 

2 

5 

0.6702 

0.3739 

0.7804 

3 

5 

3 

0.5123 

0.3535 

0.8920 

4 

1 

4 

0.6421 

0.3626 

0.8406 

4 

2 

4 

0.7262 

0.3887 

1.0014 

4 

5 

5 

0.8272 

0.3697 

1.0555 

A 

1 

1 

3.0050 

0.4416 

1.6112 

A 

2 

3 

0.5116 

2.1329 

A 

5 

1 

— 

0.4035 

1.8860 

B 

1 

2 

3.2138 

0.4057 

2.0166 

B 

2 

4 

3.4165 

0.3929 

1.8709 

B 

5 

2 

2.9061 

0.3533 

1.8516 

C 

1 

4 

3.8592 

0.7353 

3.0660 

C 

2 

4 

3.8756 

0.6058 

2.3255 

C 

5 

4 

3.6456 

0.5488 

2.5740 

D 

1 

6 

5.7705 

0.6105 

2.4467 

D 

2 

7 

4.4030 

0.5779 

2.2287 

D 

5 

4 

4.0556 

0.4792 

2.3097 
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Figure  C-9.  Optimal  Pilot  Model  Weightings  --  Roll  Axis 


Pitot  Model  Analysis 

The  dynamics  simulated  for  the  flight  test  were  evaluated  using  the  applicable 
MIL-STD-1797A  models  and  the  STI  optimal  pilot  model.  The  parameters  given  in 
Table  C-5  were  used  when  running  the  optimal  pilot  model. 

Table  C-5 

Optimal  Pilot  Model  Parameters  --  Flight  Test  Analysis 


Forcing  Function,  Yw 
(at  Aircraft  Output) 

J2 

6.25s2  +  3.54  s +1 

Motor  Noise 
Ratio,  p^ 

-20  dB 

0.08 

Visual  Indifference 
Thresholds,  T„,  Ty2 

0 

Pilot  Delay,  tp 

0.25  seconds 

Fractional  Attention 
Parameter,  f 

1 

Observation  Noise 
Ratios,  pvl  pv2 

-20  dL 

Driving  Noise 
Intensity,  Vw 

1 

The  results  of  this  analysis  are  presented  in  Tables  C-6  and  C-7.  Bode  plots  of  the 
pilot  describing  functions  predicted  by  the  STI  optimal  pilot  model  for  Cases  1  and  4  are 
shown  in  Figures  C-10  and  C-12.  Bode  plots  of  the  resulting  pilot-aircraft  systems  are 
shown  in  Figures  C-l  1  and  C-13. 
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Table  C-6 

Pitch  Axis  Pilot  Model  Predictions  --  Flight  Test  Dynamics 


Case 

C  xD2 

C-H} 

Rating 

CAP4 

SP5 

TRP6 

BW7 

Criterion 

Neal-Smith 

Criteria 

Gibson's  Criteria 

OPM8 

1 

.7  .04 

2-3 

I 

II 

■ 

II 

II 

Abrupt  Bobbling 
Tendency 

H 

2 

.4  .04 

1-2 

I 

I 

I 

I-II 

II 

Abrupt  Bobbling 
Tendency 

3.8 

3 

.7  .24 

3-5 

III 

III 

III 

III 

III 

Satisfactory 

Response 

■ 

■ 

.4  .24 

H 

III 

III 

III 

III 

III 

Satisfactory 

Response 

‘Short  Period  Damping  Ratio  5Short  Period  Criterion 

2System  Delay  'Transient  Response  Parameter 

3Cooper-Harper  Ratings  from  Flight  Test  Data  ’Bandwidth  Criterion 

4Control  Anticipation  Parameter  Criterion  “Optimal  Pilot  Model  Rating  Prediction 


Table  C-7 

Roll  Axis  Pilot  Model  Predictions  --  Flight  Test  Dynamics 


Case 

T  2 
XD 

m 

Bandwidth 

Criterion4 

Roll 

Constant 

Spiral 

Constant 

Delay 

Step 

Response 

OPM5 

D 

.4 

.04 

1-3 

2 

I 

I 

I 

I 

3.8 

B 

1 

.04 

eii 

3 

I-II 

I 

I 

I 

5.3 

C 

.4 

.24 

ISIS 

4 

I 

I 

III 

I 

5.2 

D 

1 

.24 

El 

5 

I-II 

I 

III 

I 

10 

‘Roll  Mode  Time  Constant 
JSystem  Delay 

’Cooper-Harper  Ratings  Range  from  Flight  Test  Data 
"Bandwidth  Rating  is  from  the  Regression  Formula  in  Reference  24 
’Optimal  Pilot  Model  Rating  Prediction 
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Figure  C-12.  Bode  Plot  of  Predicted  Pilot  Transfer  Function  (Case  4) 


Frequency,  f.  (radians  per  second) 

Figure  C-13.  Bode  Plot  of  Predicted  Pilot-Aircraft  System  (Case  4) 
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Appendix  D.  Frequency  Response  Data 

This  appendix  contains  the  frequency  response  analysis  plots  referred  to  in 
Chapter  4.  The  first  two  figures  display  sample  power  spectral  densities.  Figures  D-3 
through  D-10  present  representative  frequency  responses  of  stick  displacement  to  task 
error  for  each  airborne  sum-of-sines  cases.  In  these  figures,  the  top  two  plots  represent 
the  Bode  plot  of  the  pilot  and  the  bottom  plot  presents  the  coherence  as  described  in 
Chapter  4.  For  each  of  these  figures,  the  data  points  computed  by  the  frequency  response 
analysis  software  are  represented  by  asterisks.  The  dotted  lines  on  the  figures  represent 
95  percent  confidence  boundaries.  Figures  D-l  1  through  D-18  present  the 
corresponding  Bode  plots  of  the  combined  pilot-aircraft  systems.  Finally,  Figures  D-l 9 
and  D-22  are  frequency  response  plots  from  Case  1  for  the  command  path  gain  analysis 
described  in  Chapter  4. 
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Output  Power  Spectral  Density  Input  Power  Spectral  Density 


Calspan  Variable  Stability  Aircraft 
Leaijet  LJ-25,  Tail  Number  N102VS 
Date:  9- 11  Oct  93;  Pilot:  G;  Sortie  #1 
Sum-of-Sines  Tracking  Task;  Case  4;  Airborne  Data 


10' 

io‘ 

10' 

10'  . 
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Frequency,  (radians  per  second) 


Frequency,  co  (radians  per  second) 


Figure  D-l.  Power  Spectral  Density  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  4) 


D-2 
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Calspan  Variable  Stability  Aircraft 
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Figure  D-2.  Power  Spectral  Density  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  D) 
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Figure  D-3.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  1 ) 
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Figure  D-4.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  2) 
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Figure  D-5.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  3) 
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Figure  D-6.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  4) 
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Figure  D-7,  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  A) 
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Figure  D-8.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  B) 
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Figure  D-9.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  C) 
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Figure  D- 1 0.  Frequency  Response  Analysis  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  D) 
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Figure  D-l  1 .  Combined  Pilot- Aircraft  System  (Case  1 ) 
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Figure  D-12.  Combined  Pilot- Aircraft  System  (Case  2) 
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Figure  D-14.  Combined  Pilot- Aircraft  System  (Case  4) 
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Figure  D-15.  Combined  Pilot-Aircraft  System  (Case  A) 
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Figure  D-16.  Combined  Pilot-Aircraft  System  (Case  B) 
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Figure  D-17.  Combined  Pi  lot- Aircraft  System  (Case  C) 
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Figure  D-18.  Combined  Pilot-Aircraft  System  (Case  D) 
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Figure  D-19.  Gain  Investigation  Baseline  - 
Longitudinal  Stick  Deflection  to  Task  Error  (Case  1) 
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Figure  D-20.  Gain  Investigation  -  Longitudinal  Stick 
Deflection  to  Task  Error  (Case  1  with  Double  Gain) 
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Figure  D-21 .  Combined  Pilot-Aircraft  System 
--  Gain  Investigation  Baseline  (Case  1 ) 
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Figure  D-22.  Combined  Pilot-Aircraft  System 
--  Gain  Investigation  (Case  1  with  Double  Gain) 
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Appendix  E.  Numerical  Solution  Routines  and  Results 


PIFIND.M 

function  [J]-PIFTND(Ak,Bk,Ck,Ac,Bc,Cc,Ag,Bg.Cg)Q,R,Noise) 

% 

%  This  File  Finds  the  Performance  Index  Value  for  Problems 
%  With  Filter  States 
%  where 

%  Ak,Bk,Ck  -  Compensator  States 

%  Ac,Bc,Cc  *  Aircraft  States 

%  Ag,Bg,Cg  =  Disturbanc  Noise  Filter  States 

%  Q  =  State  Deviation  Weighting 

%  (Diagonal  with  #Rows-#Plant+#Filter  States) 

%  R  =  Control  Usage  Weighting 

%  (Diagonal  with  #Rows=#Rows  of  Ck) 

%  Noise  =  Noise  Matrix 

%  |Qo  0  |  Qo  =  Process  Noise  Intensity 

%  |0  Ro|  Ro 9  Measurement  Noise  Intensity 

%  (Square  with  #Cols=#Cols  of  Bg+Bk) 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

sizeak-size(Ak); 

sizeac=size(Ac); 

sizeag-size(Ag); 

sizebg-size(Bg); 

sizebk=-size(Bk); 

sizeck«size(Ck); 

A-[Ac,Cg,-Bc*Ck;zeros(sizeag(l),sizeac(l)),Ag)zeros(sizeag(l),sizeak(l));Bk*Cc,zeros(sizeak(l), 

sizeag(l)),Ak]; 

E-[zeros(sizeac(l),sizebg(2)),zeros(sizeac(l),sizebk(2));Bg,zeros(sizeag(I),sizebk(2));zeros(sizeak(l), 

sizebg(2)),Bk]; 

L=E*Noise*E'; 

X»lyap(A,L); 

m=sizeak(l); 

n-sizeac(l)+sizeag(l); 

0“sizeck(l); 
teiml=0; 
term2=0; 
for  i-l:n 

term  1  =term  1 +Q(i,i)*X(i,i); 
end 

forh«l:o 

fori=l:m 

forj-l:m 

teim2-term2+R(hdi)*Ck(h)i)*Ck(hj)*X((n+i),(n+j)); 

end;end;end 

J-terml+term2; 

return 
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PIFINDNF.M 


function  [J]=PIFINDNF(Ak,Bk,Ck,Ac,Bc,Cc,Gamma,Q,R,Noise) 

% 

%  This  File  Finds  die  Performance  Index  Value  For  Problems  with 
%  No  Filter  States 
%  where 

%  Ak,Bk,Ck  -  Compensator  States 

%  Ac,Bc,Cc  =  Aircraft  States 

%  Gamma  -  Constant  Disturbance  Noise  Matrix 

%  (Distributes  Disturbance  Noise  into  States) 

%  Q  -  State  Deviation  Weighting 

%  (Diagonal  with  #Rows=#Plant+#Filter  States) 

%  R  =  Control  Usage  Weighting 

%  (Diagonal  with  #Rows=#Rows  of  Ck) 

%  Noise  =  Noise  Matrix 

%  |Qo  0 1  Qo  =  Process  Noise  Intensity 

%  |0  Ro|  Ro  =  Measurement  Noise  Intensity 

%  (Square  with  #Cols=#Col  of  Gamma+Bk) 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

sizeak=size(Ak); 

sizeac=size(Ac); 

sizebk“size(Bk); 

sizeg“size(Gamma); 

sizeck=size(Ck); 

A=[Ac,-Bc*Ck;Bk*Cc,Ak]; 

E”[Gamma,zeros(sizeac(  l),sizebk(2));zeros(sizeak(  l),sizeg(2)),Bk] ; 

L=-E*Noise*E'; 

X=lyap(A,L); 

m-sizeak(l); 

n=sizeac(l); 

o=sizeck(I); 

term  1=0; 

term2=0; 

for  i=l:n 

term  1  =term  1  +Q(i,i)*X(i,i); 
end 

forh=l:o 
for  i=l:m 
forj=l:m 

term2=term2+R(hJh)*Ck(h,i)*Ck(hj)*X((n+i),(n+j)); 

end;end;end 

J=terml+term2; 

return 
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Figure  E-l.  Bode  Magnitude  Comparison  Plot  (First-Fourth  Order  Compensators) 
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Appendix  F.  Sub-Optimal  Pilot  Model  Algorithm  Verification 

General 

The  numerical  solution  technique  used  in  the  sub-optimal  pilot  model  is  verified 
in  this  appendix.  A  simple  example  problem  was  solved  using  both  the  standard  linear 
quadratic  Gaussian  (LQG)  approach  and  the  numerical  method  (Equations  6-19,  6-21, 
6-22,  and  6-23).  As  shown,  the  results  were  the  same  using  either  method. 


Standard  Linear  Quadratic  Gaussian  Solution 

Assume  the  controlled  element,  Yc,  and  forcing  function  filter,  Yw,  can  be 
modeled  by  the  following  transfer  functions. 


These  transfer  functions  have  the  following  state  space  representations. 


[xw\  =  [-2]  •  [xw\  +  [1]  •  u 
yw  =  [  Jl  ]  •  [xw]  +  [0]  •  u 


(F-2) 


(F-3) 


Assume  this  example  is  in  the  form  shown  in  Figure  6-2,  where  the  disturbance  occurs  at 
the  controlled  element  output  so  that  the  error  is  the  difference  between  the  filter  output 


F-l 


and  the  controlled  element  output  (e  =  Yw  -  Y^.  For  this  example,  the  desired  weighting 
on  error,  Q,  will  be  1  and  the  desired  weighting  on  control,  R,  will  be  2  in  accordance 
with  the  following  performance  index. 

GO 

J=  \[eTQe  +  uTRu]dt  (F-4) 

o 

This  example  will  assume  unit  intensity  forcing  function  driving  noise,  £(t),  and 
observation  noise,  ri(t),  such  that: 


(F-5) 

£{n(<)nr(<-T)}  =«.s(t)=8(t) 

(F-6) 

To  solve  the  standard  LQG  problem,  the  filter  states  must  be  appended  to  the 
controlled  element  states  so  that  the  system  is  in  the  following  form. 


x=Ax+Bu  +  Tl,  amd  y  =  Cx+Dy\ 

For  this  problem,  this  can  be  accomplished  using  the  following  formula: 


xc 


Ac  0 
0  Aw 


xc 

xw 


0 

Bc 


•  u  + 


0 

Bw 


y  =  [-Cc  Cw  ] 


Xc 

Xw 


+  [1]’T| 


(F-7) 


(F-8) 


where  the  optimal  LQG  control  law  is  u  =  -Kx.  Substituting  Equations  F-2  and  F-3  into 
Equation  F-8  yields  the  following  state  space  representation. 


F-2 


*C1 

’  -2 

-4 

0 

Xci 

'  1  ' 

’  0  " 

XC2 

= 

1 

0 

0 

• 

Xci 

+ 

0 

•«  + 

0 

Xw 

0 

0 

-2 

Xw 

0 

1 

y  =  [  -1  -3  Jl  ]• 


Xci 

*C2 

Xw 


+  [1]-T| 


(F-9) 


Error  can  be  weighted  using  the  following  formula: 


Q>  =  CtQC 


(F-10) 


where  Q  is  one  for  this  example,  and  the  new  state  weighting,  Q' ,  is: 


1  3  -Jl 

3  9  -3^2 

-Jl  -3^2  2 


(F- 11) 


The  noise  weighting  matrix,  Qf,  is  found  using  the  following  formula. 


Qf=r-Q0rT= 


0  0  0 
0  0  0 
0  0  1 


(F-12) 


Using  the  following  MATLAB™  command  will  return  the  optimal  LQG  compensator: 

[Ak,Bk,Ck,Dk]=lqg(A,B,C,D,W,V) 

where  A,  B,  C,  and  D  are  the  state  space  matrices  from  Equation  F-9,  W  contains  the 
weighting  matrices  in  diagonal  form  ( diag{Q R}),  and  V  contains  the  noise  matrices  in 
diagonal  form  ( diag{Qf ,  Ra }).  The  optimal  LQG  compensator  for  this  problem  is: 
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(F-13) 


-2.3570  -4.5277 

0.2670 

0 

X  = 

1 

0 

0 

■x  + 

0 

0.4313 

1.1212 

-2.5344 

0.3178 

y  =  [  0.3570  0.5277  -0.2670  ]-x 


or  in  transfer  function  form,  the  optimal  LQG  compensator,  Kj^s),  is: 

K  (s)  _  -0.0849^-0.1697^-0.3395 

LQGK  ’  s3  +  4.89 13s2  +  10.38605  +11. 1754 


(F-14) 


The  Numerical  Solution 

To  accomplish  the  numerical  solution  the  file,  PI JSOPM.M,  has  to  be  adapted 
slightly.  The  desired  compensator  form  is  now  a  second  over  a  third  order  transfer 
function  as  shown  in  the  new  file,  PIFINDER.M,  presented  at  the  end  of  this  appendix. 
Except  for  this  change,  PIFINDER.M  is  identical  to  the  routine  used  in  the  sub-optimal 
pilot  model. 

For  this  example,  the  controlled  element  state  space  matrices  remain  as  in 
Equations  F-2  and  F-3.  The  weightings,  Q  and  R,  are  the  scalar  values,  I  and  2.  Finally, 
the  Noise  matrix  is  a  diagonal  matrix  containing  the  noise  intensities,  Q0  and  R0 .  In  this 
example  that  is  a  2  by  2  identity  matrix.  Using  an  initial  compensator  coefficient  vector, 
cO,  of  all  ones,  the  following  command  will  return  the  optimal  solution. 

cmin=ftninsCpifinder',cO,[],[],Ac,Bc,Cc,Aw,Bw,Cw,Q,R,Noise) 

The  following  value  for  cmin  was  computed. 

cmin=[0.0849  0.1058  0.3218  4.0764  7.7238  10.4016] 
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Using  this  vector  as  the  new  initial  compensator  coefficient  vector,  cO,  and  running  fmins 
again  yielded: 


cmin=[0.0849  0.1697  0.3395  4.8065  10.2163  10.8360] 

This  equates  to  the  following  compensator. 

K(s)  =  0.0849,?2  +  0. 1 697s  +  0.3395 

U  j3  +  4.8065.S2  + 10.2163^+  10.8360 


(F-15) 


Notice  that  the  coefficients  are  identical,  within  the  tolerance  of  the  numerical 
search,  to  those  found  using  standard  LQG  methods.  Also  note  that  all  of  the  numerator 
coefficients  are  positive  while  those  in  Equation  F- 14  are  negative.  This  is  because  the 
standard  LQG  optimum  control  law  is  always  u  =  -Kx,  while  the  numerical  method  used 
in  the  sub-optimal  pilot  model  has  the  control  law,  u-  Kx  (see  Figure  6-2,  page  6-4). 
Thus,  the  numerical  method  used  to  minimize  the  performance  index  in  the  sub-optimal 
pilot  model  is  valid. 


MATLAB  '"File  for  Disturbance  at  Plant  Output 

function  [J]=PIFINDER(c,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise) 

% 

%  This  File  Finds  the  Minimum  Performance  Index  Value 
%  For  Problems  that  Weight  Error  through  Disturbance 
%  Noise  Injected  at  the  Aircraft  Output 
% 

%  where 

%  c  =  Compensator  Coefficients 

%  AcJBc,Cc  -  Aircraft  States 
%  Ag,Bg,Cg  =  Disturbance  Noise  Filter  States 
%  Q  =  Error  Weighting  (scalar) 

%  R  *  Control  Usage  Weighting  (scalar) 

%  Noise  =  Noise  Matrix 
%  |Qo  0 1  Qo  =  Process  Noise  Intensity 

%  |0  Ro|  Ro  =  Measurement  Noise  Intensity 

% 
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%%%%%Vo%%%%%%%%%%%%%%%%%%%Vo%%Vo%VM%%%%%%%%%%%%%%%%%%% 
%  Compensator  Form  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#/o%%%%%%%%%%%%%%%%%%%%% 
% 


% 

mim=[c(l)  c(2)  c(3)]; 
den=[l  c(4)  c(5)  c(6)]; 
den-abs(den); 

[Ak,Bk,Ck,Dk]=tf2ss(num,den); 

% 


% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Compute  J  % 

%%%%%%#/o%%%%%#/o%#/«%%%%#/«%#/o#/o%0/o%#/o%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% 

sizeak=size(Ak); 

sizeac=size(Ac); 

sizeag=size(Ag); 

A=[Ac,zeros(sizeac(l),sizeag(l)',,Bc*Ck;zeros(sizeag(l),sizeac(l)),Ag,zeros(sizeag(l),sizeak(l));-Bk*Cc 

,Bk*Cg,Ak]; 

E=[zeros(sizeac(  1  ),2);Bg,zeros(sizeag(  1 ) ,  1 )  ;zeros(sizeak(  1 ) ,  1 )  ,Bk] ; 

L=E*Noise*E'; 

X=lyap(A,L); 
m=sizeak(l); 
p=sizeac(l); 
f=sizeag(l); 
n=p^f; 
term  1=0; 
term2a=0; 
term2b=0; 
term2c=0; 
for  i=l:m 
forj=l:m 

term  1  *=term  1  +Ck(  1  ,i)*Ck(l  j)*X((n+i),(n+j)); 
end;end 
fori=l:f 
forj=l:f 

term2a=term2a+Cg(l  ,i)*Cg(l  j)*X((p+i),(p+j)); 
end;end 
fori=l:p 
forj=l:p 

term2b~term2b+Cc(l  ,i)*Cc(l  j)*X(i  j); 
end;end 
for  i=l:f 
forjml:p 

term2c=term2c+Cg(l  ,i)*Cc(l  j)*X((p+i)j); 
end;end 

J«R*term  1  +Q*term2a+Q*term2b-2*Q*term2c; 
return 
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Appendix  G.  Sub-Optimal  Pilot  Model  Data 


SOPM.M 

function  [J.Cmin, Rating]=sopm(Numc, Dene, De  lay, Numg, Deng) 

% 

%  Sub-Optimal  Pilot  Model 

% 

%  Where: 

%  Numc  =  Numerator  of  Aircraft  Dynamics 

%  Dene  -  Denominator  of  Aircraft  Dynamics 
%  Delay  -  Delay  of  Aicraft  Dynamics 
%  Numg  »  Numerator  of  Task  Forcing  Function 
%  Deng  -  Denominator  of  Task  Forcing  Function 

%  Recommend  2nd  Order  Butterworth: 

%  Numg-sqrt(2) 

%  Deng=[lA>w  sqrt(2)/bw  1] 

%  Where  'bw'  is  Butterworth  Filter  Bandwidth 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Input  Other  Variables  % 

%%%%%%%0/o0/.%%%%%%%%%#/«%%%%%#/o%%%#/«%%0/.0/o%%%0/.#/o%%%%%0/o%#/o%%%%% 

% 

if(nargin~=5),errorCIncorrect  Number  of  Input  Arguments'), end 
Q=inputCInput  Error  Weighting  (Usually  1)  =*  *); 

R-inputflnput  Control  Weighting  (Usually  4.5  Pitch,  1  Roll)  =  ■); 
rhod=input(Tnput  Observation  Noise  Ratio  (Usually  0.01)  =  ■); 

Tn=inputCMuscular  Time  Constant  (Usually  0. 1 15)  =  *); 

Pilot_Delay=inputCInput  Pilot  Delay  (sec)  (Usually  .35)  =  ■); 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Model  Parameters  are  Set  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% 

Var=|Tn, Pilot JDelay+Delay] ; 
c0=[l  1  1]; 

V-.5; 

rhodif-l; 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Performance  Index  is  Minimized  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%0/o%%%%%%%%%%%%%%%%%%%%%%% 

% 

[Ac,Bc,CcJ)c]-tf2ss(Numc,Denc); 

[Ag3g,Cg,Dg]-tf2ss<Numg,Deng); 
for  i-l:5 
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if  rhodif>0.0002  %0. 1  db  at  rtao  -  -20  dB  ( 1 0*log  1 0) 

disptflteration  ',mt2str(i)]> 

Noise-[1,0;0,V] 

Cmin=-fTninsCpLsopm',cO,Q,Q,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise,Var) 

[J(RMS]-RMS_sopm(Cmin,Ac>Bc,Cc,Ag,Bg,Cg,Q,R,Noise,Var) 

iho**V/pi/RMS 

rhodif**abs(rho-rhod); 

cO-Cmin; 

V-V*rhod/rho; 

end 

end 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Cooper-Harper  Rating  is  Computed  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

x=menuCSelect  Axis'.'Pitch'.'Roll1); 
ifx»l 

Rating--30+24 1  *log  1 0(J); 
else 

Rating«-13+117*loglO(J); 

end 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Bode  Plot  and  Transfer  Function  of  Pilot  are  Displayed  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% 

Cmin=abs(Cmin); 

K-Cmin(l); 

A”Cmin(2)/Cmin(l); 

B=Cmin(3); 

C-l/Tn; 

dispO 

dispCPilot  Parameters  are:  *) 
dispCfK  =  >um2str(K)]) 
disj^fA  =  ’,num2str(A)j) 
dispff'B  =  >>num2str(B)]) 
disp(['C  “  ’jnum2str(C)]) 
disp(['D  =  ’,num2str(Pilot_Delay)]) 

di^O 

dispCWhere  the  Pilot  is  KCs+A^exp^DsyCS+BHS+C)1) 

% 

w=-logspace(- 1 , 1 , 1 00); 

num“[Cmin(l)  Cmin(2)]; 

den=[l  Cmin(3)+1/Tn  Cmin(3)*l/Tn]; 

[mag, phase, w]=bode(num, den, w); 
mag”20*logl0(mag); 
phase-phase-w1  *  Pilot_Delay  *  1 80/pi; 
hold  off; 
axisCnormal1); 
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subplot(2U) 

semilogx(wfmag,[.  1  10], [0,0]); 
subplot(212) 

semilogx(w, phase, [.1  10],f-180 -180]) 
subplot(l  1 1) 

titlefBode  Magnitude  and  Phase  of  Pilot') 
xlabel(’Frequency  (Rad./Sec.)r) 
ylabelCPhase  (deg.)  Mag.  (dB)') 


PI SOPM.M 

function  [J)=PI_SOPM(c,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise,Var) 

% 

%  This  File  is  Called  by  SOPM.M  to  Find  the  Classical  Pilot  Model 
%  Form  that  Produces  die  Minimum  Performance  Index  Value 
% 

%  where 

%  c  -  Compensator  Coefficients 

%  Ac,Bc,Cc  -  Aircraft  States 

%  Ag,Bg,Cg  =  Disturbance  Noise  Filter  States 
%  Q  =  Error  Weighting  (scalar) 

%  R  -  Control  Usage  Weighting  (scalar) 

%  Noise  -  Noise  Matrix 

%  |Qo  0 1  Qo  -  Process  Noise  Intensity 

%  |0  Ro|  Ro  -  Measurement  Noise  Intensity 

%  VAR  -  Variables  -  [Muscular  Time  Constant,  Delay] 

% 

%%%%%%%%%%%%0/o%%%%%#/«%%%0/o%0/o#/o%%%%%0/o%%#/.0/o%0/o0/o0/o0/o%%0/o0/o0/o0/o#/o0/o%% 

%  Pilot  Model  Form  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

[num_delay,den_delay]“pade(Var(2),  1 ); 
num«conv([c(l)  c(2)],num_delay); 
den=conv([l  c(3)+l/Var(l)  c(3)*l/Var(l)],den_delay); 
den=»bs(den); 

[Ak,Bk,CkJDk]“tf2ss(num,den); 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Compute  J  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

sizeak=size(Ak); 

sizeac-size(Ac); 

sizeag=*size(Ag); 

A=[Ac^eros(sizeac(  1  ),sizeag(  1  )),Bc*Ck;zeros(sizeag(  1  ),sizeac(  1)),  Ag,zeros(sizeag(  1  ),sizeak(  1  ));-Bk*Cc 
3k*Cg,Ak]; 

E«[zeros(sizeac(l),2);Bg,zeros(sizeag(l),l);zeros(sizeak(l),l),Bk]; 

L»E*Noise*E'; 

X=lyap(A3); 
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m-sizeak(l); 
p-sizeac(l); 
ft-sizeag(l); 
n— pH-f; 
term  1=0; 
temt2a-0; 
term2b“0; 
term2c=0; 
for  i-l:m 
for  j-l:m 

term  1  -term  1  +Ck(  1  ,i)*Ck(  1  j)*X((n+i),(n+j)); 
end;end 
for  i-l:f 
for  j-l:f 

term2a**term2a+Cg(l  ,i)*Cg(l  j)*X((p+i),(p+j)); 
end;end 
for  i-l:p 
for  j“l:p 

term2b=term2b+Cc(  1  ,i)*Cc(l  j)*X(i  j); 
end;end 
for  i—1  :f 
for  j— 1  :p 

term2c-term2c+Cg(l  ,i)*Cc(l  j)*X((p+i)j); 
end;end 

J=R*term  1  +Q*term2a+Q*term2b-2*Q*term2c; 
return 


RMSJSOPM.M 

function  [J,RMS]=RMS_SOPM(c,Ac,Bc,Cc,Ag,Bg,Cg,Q,R,Noise,Var) 

% 

%  This  File  is  Called  by  SOPM.M  to  Find  the  RMS  Error 
%  and  Performance  Index  Values 
% 

%  where 

%  c  -  Compensator  Coefficients 

%  Ac3c.Cc  ”  Aircraft  States 

%  Ag,Bg»Cg  -  Disturbance  Noise  Filter  States 

%  Q  -  Error  Weighting  (scalar) 

%  R  »  Control  Usage  Weighting  (scalar) 

%  Noise  -  Noise  Matrix 

%  (Qo  0 1  Qo  -  Process  Noise  Intensity 

%  |0  Ro)  Ro  -  Measurement  Noise  Intensity 

%  VAR  -  Variables  -  [Muscular  Time  Constant,  Delay] 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Pilot  Model  Form  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

[num_delay,den_delay]=pade(Var(2),  1); 
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num«conv([c(l)  c(2)],num_delay); 

den-conv([l  c(3)+l/Var(l)  c(3)*l/Var(l)],den_delay); 

den-abs(den); 

[Ak,Bk,Ck,Dk]=tf2ss(num,den); 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Compute  J  % 

%%%%%%%%%%%%%°/o%%%%%%%%%%0/ i%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

sizeak=size(Ak); 

sizeac”=size(Ac); 

sizeag-size(Ag); 

A=[Ac,zeros(sizeac(  l),sizeag(l)),Bc*Ck;zeros(sizeag(  l),sizeac(  1)),  Ag,zeros(sizeag(  1  ),sizeak(  l));-Bk*Cc 
,Bk*Cg,Ak]; 

E“[zeros(sizeac(  1  ),2);Bg,zeros(sizeag(  1 ),  1  );zeros(sizeak(  1 ),  1  ),Bk] ; 

L=E*Noise*E'; 

X=lyap(A,L); 

m=sizeak(l); 
p-sizeac(l); 
f-sizeag(l); 
n=p+f; 
terml“0; 
term2a=0; 
term2b-0; 
term2c=0; 
for  i*l:m 
for  j=l:m 

terml=terml+Ck(l  ,i)*Ck(l  j)*X((n+i),(n+j)); 
end;end 
fori-l:f 
forj=l:f 

term2a=term2a+Cg(  1  ,i)*Cg(l  j)*X((p+i),(p+j)); 
end;end 
for  i*l:p 
for  j—l:p 

term2b=term2b+Cc(l  ,i)*Cc(  1  j)*X(i  j); 
end;end 
fori-l:f 
for  j-l:p 

term2c=term2c+Cg(  1  ,i)*Cc(l  j)*X((p+i)j); 
endjend 

J=R*term  1  +Q*term2a+Q*term2b-2*Q*term2c; 

RMS*term2a+term2b-2*term2c; 

return 
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Phase  (degrees) 


Figure  G-l.  Bode  Magnitude  Plot  of  Predicted  Pilot  Describing  Functions 
--  Sub-Optimal  Pilot  Model  Cases  1  Through  4 


Figure  G-2.  Bode  Phase  Plot  of  Predicted  Pilot  Describing  Functions 
—  Sub-Optimal  Pilot  Model  Cases  1  Through  4 


Figure  G-3.  Bode  Magnitude  Plot  of  Predicted  Pilot  Describing  Functions 
—  Sub-Optimal  Pilot  Model  Cases  A  Through  D 


Figure  G-4.  Bode  Phase  Plot  of  Predicted  Pilot  Describing  Functions 
—  Sub-Optimal  Pilot  Model  Cases  A  Through  D 
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